AFRL-VS-HA-TR-98-0086 


Validation  of  Ionospheric  Models 


Dwight  T.  Decker 
Patricia  H.  Doherty 


Boston  College 

Institute  for  Scientific  Research 
Chestnut  Hill,  MA  02467 


30  April  1997 


Scientific  Report  No.  1 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  IS  UNLIMITED. 


AIR  FORCE  RESEARCH  LABORATORY 
Space  Vehicles  Directorate 
29  Randolph  Rd 

AIR  FORCE  MATERIEL  COMMAND 
Hanscom  AFB,  MA  01731-3010 


20030902  060 


This  technical  report  has  been  reviewed  and  is  approved  for  publication. 


*  (This  report  has  been  reviewed  by  the  ESC  Public  Affairs  Office  (ESC/PAM)  and  is  releasable 
to  the  National  Technical  Information  Service  (NTIS).) 


Qualified  requestors  may  obtain  additional  copies  from  the  Defense  Technical  Information 
Center  (DTIC).  *  (All  others  should  apply  to  the  NTIS.) 


If  your  address  has  changed,  or  if  you  wish  to  be  removed  from  the  mailing  list,  or  if  the 
addressee  is  no  longer  employed  by  you,  please  notify  AFRL/VSIM,  29  Randolph  Road, 
Hanscom  AFB,  MA  01731-3010.  This  will  help  us  to  maintain  a  current  mailing  list. 


DO  NOT  RETURN  COPIES  OF  THIS  REPORT  unless  contractual  obligations  or  notices  on  a 
specific  document  require  that  it  be  returned. 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


1 


Public  reoorting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this 
collection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services.  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson 
Davis  Highway.  Suite  1204,  Arlington,  VA  22202-4302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project  (0704-0188),  Washington,  DC  20503. 

1.  AGENCY  USE  ONLY  (Leave  blank)  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

30  April  1997  Scientific  Report  No.  1 

4.  TITLE  AND  SUBTITLE 

VALIDATION  OF  IONOSPHERIC  MODELS 

5.  FUNDING  NUMBERS 

Contract: 

F19628-96-C-0039 

PE61102F 

PR  1010 

TAIM 

WU  AC 

6.  AUTHOR(S) 

Dwight  T.  Decker 

Patricia  H.  Doherty 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Boston  College 

Institute  for  Scientific  Research 

140  Commonwealth  Avenue 

Chestnut  Hill,  MA  02167-3862 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING /MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Research  Laboratory 

29  Randolph  Road 

Hanscom  AFB,  MA  01731-3010 

Contract  Manager:  John  Retterer/VSBP 

10.  SPONSORING /MONITORING 

AGENCY  REPORT  NUMBER 

AFRL-VS-HA-TR-98-0086 

11.  SUPPLEMENTARY  NOTES 

12a.  DISTRIBUTION /AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

12b.  DISTRIBUTION  CODE 

13.  ABSTRACT  (Maximum  200  words) 

:  1  * 

In  the  last  year,  we  have  studied  several  issues  involving  the  validation  of  ionospheric  models 

and  the  use  of  observational  databases.  Modeling  work  has  included  a  mid-latitude 

1  intereomparison  of  five  first  principal  physical  models  for  the  earth's  ionosphere.  Work  at  low 
latitudes  has  involved  using  the  Parameterized  Ionospheric  Model  (PIM)  to  improve  the 
International  Reference  Ionosphere  (IRI90)  low  latitude  ionospheric  specification.  At  high 
latitudes,  the  work  has  involved  attempts  to  produce  blobs  at  the  times  and  locations  of  actual 
observations.  The  usefulness  of  GPS  occultations  have  been  studied  using  a  series  of  simulations 
!  of  such  occultations  and  examining  inversions  of  the  simulations.  Other  studies  have  been  made 
on  various  aspects  of  using  GPS  measurements  for  ionospheric  research.  These  include  the  solar 
cycle  dependence  of  the  Total  Electron  Content  (TEC),  the  effects  of  anti-spoofing,  the  TEC  of 
the  protonsphere,  the  effects  of  instrumental  bias,  and  spatial  and  temporal  variations  in  the  TEC. 

14.  SUBJECT  TERMS  — 

Ionospheric  weather,  Global  F-region  modeling,  Total  electron  content^ 
(TEC),  GPS,  GPS/MET,  Occultations,  Plasma  structure,  Blobs 

15.  NUMBER  OF  PAGES 

16.  PRICE  CODE  ; 

17.  SECURITY  CLASSIFICATION 

OF  REPORT 

UNCLASSIFIED 

18.  SECURITY  CLASSIFICATION 

OF  THIS  PAGE 

UNCLASSIFIED 

19.  SECURITY  CLASSIFICATION 

OF  ABSTRACT 

UNCLASSIFIED 

20.  LIMITATION  OF  ABSTRACT 

SAR 

NSN  7540-01-280-5500  Standard  Form  293  (Rev.  2-89) 


Pr^sr-ibed  by  ANSI  Std  739-18 


,’98-1)2 


TABLE  OF  CONTENTS 


Page 

1.  INTRODUCTION .  1 

2.  MODELING  STUDIES .  1 

2.1  PRIMO .  1 

2.2  Boundary  and  Subauroral  Blobs .  2 

2.3  Simulations  of  GPS/MET  Ionospheric  Observations .  2 

2.4  Improving  IRI90  Low-Latitude  Ionospheric  Specification .  19 

2.5  Ionospheric  Modeling  Using  Systems  of  Plasma  Kinetic  Equations -  37 

3.  OBSERVATIONAL  DATABASES .  38 

3.1  Solar  Cycle  Dependence  of  Absolute  Ionospheric  Total  Electron 

Content .  38 

3.2  Absolute  Real-Time  Ionospheric  Measurements  from  GPS 

Satellites  in  the  Presence  of  Anti-Spoofing .  39 

3.3  Total  Electron  Content  of  the  Earth's  Protonosphere .  39 

3.4  Instrumental  Bias  in  GPS  Measurements .  41 

3.5  Statistics  of  the  Spatial  and  Temporal  Variations  in  Ionospheric 

Range  Delay .  48 

3.6  AirglowData .  49 

3.7  Signature  of  Equatorial  Scintillation  Onset  in  Near-Sunset  DMSP 

Satellite  Data .  49 

4.  REFERENCES .  51 

5.  PRESENTATIONS  AND  PROCEEDINGS .  54 

6.  JOURNAL  ARTICLES  .  64 

Modeling  the  Formation  of  Polar  Cap  Patches  Using  Large  Plasma 

Flows . 65 

iii 


Collisional  Degradation  of  the  Proton-H  Atom  Flues  in  the  Atmosphere: 

A  Comparison  of  Theoretical  Techniques .  gg 

Total  Electron  Content  Over  the  Pan-American  Longitudes: 

March-April  1994  .  jqq 


iv 


1.  INTRODUCTION 


The  objective  of  this  research  is  to  obtain  various  ionospheric  measurements  from  a  wide  range 
of  geographic  locations  and  to  utilize  the  resulting  databases  to  validate  the  theoretical 
ionospheric  models  that  are  the  basis  of  the  Parameterized  Real-time  Ionospheric  Specification 
Model  (PRISM)  and  the  Ionospheric  Forecast  Model  (IFM).  Thus  our  various  efforts  can  be 
generally  categorized  as  either  modeling  studies  or  work  concerning  observational  databases. 


2.  MODELING  STUDIES 

2.1.  PARAMETERIZED  REAL-TIME  IONOSPHERIC  MODELS  AND 
OBSERVATIONS  WORKSHOP  (PRIMO) 

Over  the  last  several  years,  we  have  participated  in  the  PRIMO  workshop  at  the  annual  CEDAR 
meeting.  This  PRIMO  working  group  has  focused  on  the  intercomparison  of  five  first  principles 
physical  ionospheric  models  and  a  manuscript  has  been  written  describing  the  results.  While 
overall  the  various  models  compared  well,  there  were  some  remaining  differences  that  were 
disturbing.  In  particular,  there  were  differences  between  the  Global  Theoretical  Ionospheric 
Model  (GTM)  and  the  Field  Line  Interhemi  spheric  Plasma  (FLIP)  model  that  appeared 
inconsistent. 

As  a  result,  we  repeated  the  four  runs  that  were  originally  done  for  the  PRIMO  paper.  For  all 
four  runs,  we  used  the  HWM90  wind  model  rather  than  the  HWM87  model  to  be  consistent  with 
what  was  used  by  FLIP  and  Time  Dependent  Ionospheric  Model  (TDIM).  We  further  discovered 
that  the  original  two  solar  maximum  runs  had  not  been  done  for  precisely  the  same  days  as  stated 
in  the  paper.  This  had  little  effect  on  the  December  case;  however,  for  the  June  case,  it  made  a 
difference  in  the  magnetic  conditions  that  led  to  changes  as  significant  as  those  caused  by 
changing  wind  models.  We  also  decided  to  reduce  our  nocturnal  source  of  ionization  to  its 
“physical  value”.  In  previous  runs,  we  had  used  an  enhanced  nocturnal  source  to  ensure  that  the 
nighttime  densities  during  winter  solstice  were  maintained  at  a  realistic  level.  However,  for  the 
purposes  of  model  comparisons  this  simply  confused  the  issue  and  obscured  under  what 
conditions  GTIM  has  problems  maintaining  the  nighttime  ionosphere. 

These  various  changes  have  produced  a  more  consistent  comparison  between  GTIM  and  FLIP. 
For  the  winter  cases,  GUM's  nmf2  goes  from  being  a  bit  below  FLIP  at  night  to  a  bit  above  in  the 
daytime;  while  GTIM's  h„f2  is  generally  10  -15  km  above  FLIP.  The  summer  cases  have  GTIM 
h^2  going  from  around  5  km  below  FLIP  at  night  to  5-10  km  above  during  the  daytime.  GTIM's 
nnf2  generally  runs  a  bit  above  FLIP  for  the  summer  runs. 

We  also  performed  some  runs  with  the  neutral  wind  reduced  by  50  percent  as  was  done  in  the 
TDIM  runs.  This  improved  some  of  the  comparisons  between  GTIM  and  TDIM;  however,  in 
some  cases  it  made  the  comparisons  worse.  Discussions  with  Jan  Sjoka  indicated  that  they  also 
noted  something  odd  about  some  of  the  TDIM  runs.  Accordingly,  we  have  planned  no  more 
efforts  on  this  issue  until  they  understand  the  details  of  what  is  going  on  with  these  particular 
TDIM  runs. 
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2.2.  Boundary  and  Subauroral  Blobs 


In  previous  work,  we  have  shown  that  GTIM  can  simulate  F-region  electron  density 
enhancements  known  as  boundary  blobs.  In  that  work  we  did  not  compare  to  any  particular  data, 
rather  we  produced  structures  that  simply  had  many  of  the  features  associated  with  boundary 
blobs.  The  next  step  was  to  modify  the  simulation  to  see  if  blobs  could  be  produced  at  the  times 
that  they  have  been  observed.  It  was  not  obvious  that  this  could  be  done  because  our  first 
simulations  would  not  produce  blobs  at  the  time  appropriate  for  viewing  by  the  Chatanika  radar. 
We  chose  to  simulate  the  conditions  of  an  experiment  carried  out  on  29  January  1979  during 
which  the  Chatanika  radar  observed  blobs  for  over  six  hours  in  the  evening  sector.  Our  previous 
simulation  was  modified  in  two  ways.  First,  our  variation  of  the  convection  pattern  was  shifted  in 
UT  so  as  to  have  the  potential  of  producing  blobs  at  the  appropriate  time.  Second  we  shifted 
from  December  solstice  to  29  January.  Both  of  these  changes  were  crucial.  The  shift  away  from 
solstice  was  crucial  for  allowing  high  density  plasma  to  be  transported  from  the  daytime  across 
the  polar  cap  into  the  auroral  regions  at  the  appropriate  time.  The  shift  in  time  of  when  the 
convection  pattern  was  varied  ensured  that  this  plasma  was  directed  at  the  appropriate  time  into 
the  evening  sector.  The  resultant  blobs  appeared  at  the  appropriate  time  and  were  comparable  in 
magnitude  to  what  was  actually  seen  in  the  experiment.  The  next  step  will  be  to  change  the  size 

of  the  convection  pattern  to  see  if  we  can  produce  the  blobs  both  at  the  correct  time  and  correct 
latitude. 

2.3.  Simulations  of  GPS/MET  Ionospheric  Observations 

GPS/MET  is  a  program  using  a  low-earth-orbiting  (LEO)  satellite  instrumented  with  a  GPS 
receiver  to  perform  limb  soundings  of  the  atmosphere  using  radio  occultation  of  GPS  satellites. 
The  primary  purpose  of  the  program  is  to  perform  remote  sensing  of  the  troposphere  and 
stratosphere.  However,  the  limb-viewing  geometry  does  provide  a  unique  opportunity  for 
ionospheric  research.  The  useful  measured  quantity  in  this  case  is  the  Total  Electron  Content 
(TEC)  between  the  low-earth-orbiting  satellite  and  a  particular  GPS  satellite.  In  contrast  to 
ground-based  GPS  measurements,  the  limb  viewing  geometry  can  provide  the  type  of  vertical 
information  that  is  poorly  measured  by  ground-based  GPS  measurements. 

The  issue  of  interest  is  that  of  how  well  we  can  invert  such  measurements  of  integrated  electron 
density  into  the  electron  density  itself  as  a  function  of  1,  2,  or  3  dimensions.  Previous  workers 
have  shown  that  for  the  case  of  two  dimensions  a  single  occultation  does  not  uniquely  define  the 
electron  density.  Rather,  one  has  to  deal  with  an  underdetermined  system  of  equations.  One 
approach  to  this  problem  has  been  to  apply  a  two  dimensional  tomographic  formalism  to 
inverting  a  single  occultation  to  produce  a  single  electron  density  profile  (EDP)  by  assuming  a 
particular  horizontal  ionospheric  behavior.  This  idea  has  been  applied  to  cases  assuming 

spherical  symmetry  and  cases  with  a  horizontal  gradient  that  is  the  same  at  all  altitudes  \Haii  et 
al.,  1994].  JJ 

We  are  working  on  assessing  how  useful  such  an  approach  can  be.  Our  first  goal  is  to  quantify 
how  well  an  inversion  can  work  under  various  ionospheric  conditions  when  we  assume  spherical 
symmetry.  We  have  begun  by  writing  an  implementation  of  the  Abel  inversion.  In  Figure  1,  we 
show  the  geometry  of  an  ideal  occultation.  A  low-earth-orbiting  satellite  receives  a  GPS  signal 
which  can  be  occulted  by  the  earth.  The  occultation  is  assumed  to  occur  in  a  plane  and  the 
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FIGURE  1.  GPSMET  Ionospheric  Occulation 


motion  of  the  GPS  satellite  is  neglected.  By  subtracting  the  signals  taken  along  the  same  ray 
path,  the  TEC  below  the  altitude  of  the  LEO  satellite  can  be  determined.  In  order  to  test  the 
usefulness  of  the  Abel  inversion,  we  have  used  an  ionospheric  model  to  simulate  the  TEC  that 
would  be  measured  in  such  a  geometry,  inverted  the  simulated  observation,  and  then  compared  to 
the  profile  given  by  the  ionospheric  model.  The  model  used  has  been  the  Parameterized 
Ionospheric  Model  (PIM).  In  Figure  2,  we  give  an  example  for  a  solar  maximum  equinox  of  a 
simulated  GPS/MET  TEC  observation  using  a  spherically  symmetric  ionosphere  based  on  the 
PIM  during  magnetically  quiet  conditions.  The  center  of  the  occultation  region  is  at  0° 
geographic  latitude  and  longitude  and  is  for  a  local  time  of  midnight.  Figure  3  shows  the  profile 
derived  from  the  Abel  inversion  (short  dashes)  and  the  original  profile  from  the  PIM  (long 
dashes).  This  simply  shows  how  well  the  inversion  can  work  when  you  have  a  spherically 
symmetric  ionosphere. 

The  next  step  has  been  to  perform  a  series  of  simulations  using  the  full  PIM  ionosphere,  that  is, 
an  ionosphere  that  has  not  been  constructed  to  be  symmetric.  In  Figure  4,  we  show  a  daytime 
case  with  the  center  of  the  occultation  at  -40°  North  and  0°  longitude.  Again,  the  short  dashes 
are  the  Abel  inversion  and  the  long  dashes  are  the  PIM  profile.  We  see  that  in  the  daytime  mid¬ 
latitudes  the  Abel  can  work  very  well.  Plotted  on  the  right  is  the  fractional  difference  of  the  two 
profiles.  The  .1  and  .3  levels  are  marked  by  vertical  lines  about  zero.  We  see  that  the  agreement 
is  better  than  10  percent  above  200  km.  Figure  5  shows  the  same  comparisons  except  for  10°N 
latitude.  Again,  in  the  F-region,  the  two  profiles  agree  to  within  10  percent.  However,  the  E 
region  shows  large  differences  with  the  Abel  inversion  producing  negative  densities.  The  case 
for  40°  North  is  shown  in  Figure  6,  where  we  see  good  agreement  in  the  E  region  but  the  F 
region  agreement  is  much  poorer.  We  summarize  our  daytime  results  in  Figure  7  where  we  show 
contours  of  fractional  differences  versus  altitude  and  latitude.  In  general,  the  F  regions  agree  to 
better  than  20  percent.  The  one  exception  is  around  40°  North  where  there  is  a  strong  F  region 
gradient  in  the  PIM  ionosphere  that  causes  the  disagreement.  The  E  region  is  clearly  a  problem 
with  the  greatest  difficulties  occurring  around  the  geomagnetic  equator  (10°  North  geographic). 

Turning  to  the  nighttime,  we  show  in  Figure  8  results  for  2000  local  time  at  -30°  North.  We  see 
fair  agreement  in  the  F-region  (30  percent)  and  very  poor  agreement  in  the  E  region.  Figure  9 
gives  the  10°  North  results  and  shows  poor  agreement  at  all  altitudes.  In  Figure  10  we  again 
summarize  by  showing  contours  of  fractional  differences  versus  latitude  and  altitude.  We  see 
that  at  night  the  Abel  inversion  does  a  much  poorer  job  than  in  the  daytime. 

To  this  point  all  the  simulated  occultations  have  been  done  assuming  a  north-south  viewing 
direction.  We  have  also  tried  some  east-west  viewing  simulations  for  the  daytime  ionosphere. 
Figures  11-13  show  three  cases  for  -60°,  0°,  and  20°  North  respectively.  Figure  14  shows 
contours  of  fractional  differences  versus  latitude  and  altitude.  We  see  that  in  this  case  the  Abel 
does  an  excellent  job  throughout  the  F  region  (above  200  km).  Future  work  on  this  topic  will 
consist  of  performing  these  simulations  for  a  variety  of  geophysical  conditions. 


The  Abel  inversion  described  above  is  an  example  of  using  what  is  known  as  continuous  inverse 
theory.  The  main  limitation  of  the  Abel  is  that  it  assumes  a  very  particular  behavior  for  the 
horizontal  behavior  of  the  ionosphere.  We  have  developed  an  alternative  approach  that  uses 
discrete  inverse  theoiy  for  performing  the  inversion.  In  this  approach,  the  integral  equation  that 
defines  the  relationship  between  electron  density  and  the  TEC  measurements  is  converted  into  a 
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FIGURE  3.  PIM  and  Abel  Inverted  Profile  Sph.  Sym.  PIM 
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PIM  and  Abel  profiles  PIM  -40  +0  +0  12ut  PIM  and  Abel  differences  PIM  -40  +0  +0  12ut 


fractional  difference 


and  Abel  differences  PIM  +10+0+0  12ut 


FIGURE  5.  (a)  PIM  and  Abel  Profi  les  PIM  + 1 0  _0  _0  1 2  UT 

(b)  PIM  and  Abel  Differences  PIM  +10+0+0  12  UT 
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FIGURE  7.  Fractional  Difference  0  Long  0  Azimuth  12  UT 
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FIGURE  9.  (a)  PIM  and  Abel  Profiles  PIM  +10  +0  +0  20  UT 

(b)  PIM  and  Abel  Differences  PIM  +10  +0  +0  20  UT 
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FIGURE  10.  Fractional  Difference  0  Long  0  Azimuth  20  UT 
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FIGURE  14.  Fractional  Difference  0  Long  0  Azimuth  14  UT 
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matrix  equation,  which  is  then  solved  using  various  techniques  of  linear  algebra  [Hajj  et  al., 
1994],  The  main  advantage  of  this  matrix  equation  solution  (MES)  approach  is  that  while  we 
may  not  know  the  horizontal  behavior  of  the  ionosphere  it  does  give  a  method  that  can  treat 
whatever  horizontal  dependence  we  may  choose.  It  is  also  the  approach  that  has  the  flexibility  to 
handle  a  priori  constraints  that  are  used  in  tomography  to  reduce  the  effects  of  what  are  often  ill- 
posed  problems. 


In  developing  a  MES  method,  some  time  was  spent  assessing  whether  any  horizontal  information 
could  be  inferred  from  the  TEC  measurements  by  accepting  a  poorer  description  of  the  altitude 
dependence  of  the  EDP.  However,  the  underdetermined  nature  of  the  problem  produced  unstable 
and  incorrect  solutions.  We  did  find  that  the  MES  method  can  give  similar  results  to  the  Abel 
when  it  is  assumed  that  the  horizontal  dependence  is  a  constant.  Further,  if  the  actual  horizontal 
dependence  of  the  ionosphere  is  assumed  then  the  actual  EDP  at  the  center  of  the  occultation 
region  can  be  recovered.  Our  implementation  of  the  method  is  not  yet  as  robust  as  desired 
because  in  the  case  of  nighttime  profiles,  where  the  E  region  can  be  three  orders  of  magnitude 
smaller  than  the  F2  peak  density,  the  method  has  trouble  reproducing  the  E  region  even  when  the 
actual  horizontal  behavior  is  provided.  We  believe  this  arises  from  discretization  errors  and  that 
improvements  to  the  numerics  of  our  algorithm  may  be  able  to  cure  this  problem. 

We  have  used  the  MES  approach  to  explore  the  possibility  of  using  other  measurements  to  help 
constrain  the  inversion  process.  In  Figure  15,  we  show  two  examples  of  using  an  in-situ  electron 
density  measurement  to  help  constrain  the  inversion.  In  the  left  panel,  we  show  a  PIM  profile 
(long  dashes),  an  Abel  inversion  of  simulated  GPS/MET  observations  (short  dashes),  and  a  MES 
inversion  (dots)  that  has  applied  to  all  altitudes  the  horizontal  dependence  of  the  ionosphere  that 
would  be  observed  by  an  in-situ  satellite  measurement.  In  the  right  panel,  we  show  the  same 
PIM  and  Abel  profiles  but  the  third  profile  is  a  MES  inversion  using  the  in-situ  observed 
horizontal  dependence  down  to  the  Abel  determined  /i,,,/)  and  a  constant  horizontal  behavior  at 

lower  altitudes.  The  main  effect  of  having  the  in-situ  measurement  was  to  improve  the  quality  of 
agreement  in  the  topside  of  the  EDPs. 


Work  has  continued  on  the  Abel  transform  approach,  we  have  started  with  to  assess  the  impact  of 
measurement  errors  on  the  transform.  In  Figure  16,  we  show  in  the  left  panel  two  EDPs  and  in 
the  right  panel  the  fractional  differences  between  the  two  profiles.  The  curve  of  long  dashes  in 
the  left  panel  is  a  PIM  profile  from  a  symmetric  ionosphere.  The  curve  of  short  dashes  is  the 
Abel  inverted  profile  when  a  3  TEC  unit  error  is  added  to  the  simulated  GPS/MET  observations. 
We  see  that  the  F-region  profile  that  is  recovered  is  fairly  insensitive  to  such  errors  whereas  the 
small  E-region  densities  are  quite  sensitive  to  them.  This  suggests  that  the  inversion  technique 
may  be  quite  robust  for  producing  good  F  region  and  topside  measurements  but  as  a  technique 
for  nighttime  E-region  it  may  be  quite  problematic. 

We  have  also  made  our  first  two  attempts  at  inverting  actual  GPS/MET  observations.  We  chose 
two  occultations  that  occurred  on  Day  298,  1995.  The  first  case  involved  GPS  #16  at  around 
13:49  UT/10:15  LT  and  the  occultation  region  was  near  Sao  Paulo,  Brazil.  In  Figure  17,  we 
compare  the  actual  TEC  profile  (solid  curve)  and  the  PIM  simulation  (dotted  curve)  for  this  case. 
Given  that  we  are  comparing  a  particular  day's  measurement  to  the  climatology  in  the  PIM,  the 
differences  are  not  unreasonable.  In  Figure  18,  we  show  in  the  left  panel  the  Abel  inversion  of 
the  data  (solid  curve)  and  the  PIM  EDP  (dashed  curve)  for  these  conditions.  In  the  right  panel, 
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we  again  show  the  Abel  inversion  but  compared  to  a  Bent  Model  EDP.  The  differences  seen 
between  the  Abel  derived  F2  peak  and  the  PIM  or  Bent  F2  peaks  is  not  unreasonable  considering 
the  day-to-day  variability  of  the  ExB  vertical  drift  that  is  observed  at  low  latitudes.  The  second 
case  involved  GPS  #6  at  around  21:23  UT/  22:44  LT  and  the  occultation  region  was  around  San 
Vito,  Italy.  Figure  19  shows  the  observed  TEC  profile  (solid  curve)  and  the  PIM  (dotted  curve) 
simulation  for  this  case.  Again  the  simulated  TEC  is  reasonable  compared  to  the  observed.  In 
Figure  20,  we  show  in  the  left  panel  the  Abel  inversion  of  the  data  (solid  curve)  and  the  PIM 
EDP  (dashed  curve).  The  right  panel  contrasts  the  Abel  inversion  (solid  curve)  and  the  Bent 
Model  EDP  (dashed  curve).  The  F2  peaks  again  are  in  reasonable  agreement;  however,  the  PIM 
EDP  topside  is  very  different  from  the  Abel  inferred  profile.  Our  first  thought  was  that  this 
might  be  an  indication  of  some  horizontal  structure  that  is  ignored  in  the  Abel  inversion.  On  the 
other  hand,  the  PIM  topside  consists  only  of  0+’  and  DMSP  observations  indicate  that  at  this  time 
and  location  Ff  may  very  well  be  the  dominant  species  at  around  800  km.  The  Bent  Model  on 
the  other  hand  is  an  empirical  model  that  should  reflect  any  observable  solar  minimum  FT 
morphology.  In  the  right  panel,  we  see  that  the  Bent  profile  does  show  signs  of  FT  dominance 
above  600  km  and  this  suggests  that  the  shape  of  the  Abel  inferred  profile  is  reflecting  a 
transition  of  0+  to  HP  dominance  at  around  500  km.  This  raises  the  possibility  that  the 
GPS/MET  database  could  become  an  excellent  source  of  data  for  studying  the  0+  to  FP  transition 
height  during  solar  minimum  conditions. 

2.4.  Improving  IRI90  Low-Latitude  Ionospheric  Specification 

In  the  middle  of  the  year,  we  submitted  a  paper  to  Radio  Science  on  this  topic.  Two  referee 
reports  were  sent  to  us  and  considerable  work  was  done  to  respond  to  those  reports.  We  began 
by  performing  a  more  extensive  review  of  the  literature  led  to  provide  a  more  extensive 
discussion  of  IRI,  its  successes,  and  its  failures.  This  review  leads  to  our  writing  the  following 
extended  introduction  to  the  paper. 

Introduction 


The  International  Reference  Ionosphere  (IRI)  is  a  global  empirical  model  that  specifies  the 
monthly  average  of  the  electron  density,  electron  temperature,  ion  temperature,  and  ion 
composition  from  80  km  to  1000  km.  It  has  been  developed  as  a  joint  URSI/COSPAR  project 
and  first  appeared  as  tables  of  profiles  presented  at  the  XVII  General  Assembly  of  URSI  in  1972. 
However,  the  general  philosophy  was  to  develop  a  computer  model  based  on  critically  evaluated 
observations  and  so  IRI  quickly  evolved  from  a  set  of  tables  to  a  computer  program  (IRI-78) 
available  to  the  research  community.  Over  the  years,  testing  and  modification  of  IRI  has 
continued  with  extensive  participation  by  the  international  research  community.  The  result  is 
that  IRI  has  evolved  and  improved  through  several  versions  (IRI-80,  IRI-86,  IRI-90,  IRI-95)  and 
has  become  a  valuable  research  tool  [Bilitza  et  al.,  1993;  Szuszczewicz  et  al.,  1996]. 

As  a  data-based  model,  the  IRI  is  naturally  only  as  good  as  the  available  observations.  In  the 
case  of  electron  density,  the  data  sources  are  ionosonde  measurements,  incoherent  scatter 
measurements,  rocket  measurements,  and  a  selection  of  satellite  measurements  [Bilitza,  1990]. 
Because  the  F2-peak  parameters  (peak  density  and  height)  are  based  on  ground  based  ionosonde 
measurements,  which  are  for  the  most  part  located  in  the  mid-latitudes,  there  is  a  lack  of 
observations  from  the  low  latitude  regions,  the  high  latitudes/polar  cap  and  the  ocean  areas. 
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Likewise,  the  thickness  of  the  bottomside  F-region  is  based  on  ionosonde  measurements  and  has 
similar  problems  in  coverage.  The  topside  F-region  is  based  on  Alouette  1  and  2  topside  sounder 
measurements  from  1962  -  1966  [Llewellyn  and  Bent,  1973],  however  this  period  included  only 
solar  minimum  and  low  solar  moderate  (F10.7  <  130)  conditions.  Given  this  non-uniform 
coverage  in  both  location  and  time  of  the  various  databases,  it  is  not  surprising  that  IRI  does  an 
excellent  job  of  specifying  the  ionosphere's  climatology  for  some  regions  and  times  while  havin° 
problems  with  other  regions  and  times. 

In  this  paper,  we  focus  on  the  low  latitude  F2-region,  a  region  where  IRI  has  had  some  success 
reproducing  the  general  climatology  especially  for  lower  levels  of  solar  flux,  but  where  there  are 
some  persistent  disagreements  between  model  and  observations  [Abdu  et  al„  1996],  Early  work 
based  on  total  electron  content  (TEC)  measurements  found  examples  at  low  latitudes  where  IRI 
"severely  underestimated"  the  observed  TEC  "during  both  day  and  night."  Analysis  of  these 
cases  suggested  that  the  problem  was  in  the  relative  shape  of  IRI's  F2-region  at  low  latitudes 
[Bilitza,  1985].  The  topside  of  this  layer  comes  from  the  Bent  model  [Llewellyn  and  Bent,  1973] 
and  a  correction  to  IRI's  implementation  of  that  model  appeared  in  IRI-86.  While  this  produced 
better  agreement  with  satellite  and  incoherent  scatter  data,  it  was  noted  at  the  time  that  the  new 
topside  formula  "only  removed  half  of  the  difference  between  measured  and  predicted"  TEC  at 
the  magnetic  equator  [Bilitza,  1985],  Further,  while  IRI-86  appeared  to  do  a  reasonable  job  at 
low  solar  activity,  comparisons  with  in  situ  satellite  measurements  during  high  solar  activity 
indicated  that  IRI's  topside  profile  still  significantly  underestimated  the  observations  [Bilitza  et 
al,  1987],  The  other  half  of  the  "shape  of  the  F2-region"  is  the  bottomside  of  the  F2-layer,  which 
is  described  by  an  analytic  function  parameterized  in  terms  of  a  thickness  parameter  for  which 
IRI  offers  two  options.  The  recommended  option  in  most  cases  comes  from  the  Gulyaeva  et  al. 
[1987]  model  for  the  half-density  height.  However,  this  bottomside  thickness  has  also  had  some 
problems  at  low  latitudes.  Mahajan  et  al.  [1995]  showed  using  Arecibo  incoherent  scatter  radar 
data  that  IRI  produced  a  bottomside  that  was  too  thick  by  differing  amounts  depending  on 
season  and  improvements  to  IRI  have  been  suggested  based  on  this  work  [Gulyaeva  et  al.,  1996], 
Comparisons  to  digisonde  data  [Reinisch  and  Huang,  1996]  have  likewise  shown  some 
differences  between  data  and  IRI  with  the  most  dramatic  differences  occurring  in  August  1993  at 
Jicamarca  where  IRI's  bottomside  F2-layer  was  seen  to  be  tens  of  kilometers  thinner  than 
observed.  Finally,  we  turn  to  the  F2-peak  density  (NmF2)  and  height  ( hmF2 )  which  are  used  in  IRI 
to  complete  the  definition  of  the  F-region.  For  these  two  parameters,  the  most  persistent  problem 
at  low  latitudes  has  been  with  hmF2  and,  in  particular,  the  inability  of  IRI  in  certain  periods  to 
reproduce  the  observed  rising  in  altitude  of  the  F2-layer  just  after  sunset  [Bilitza  et  al.  1993- 
Abdu  et  al.,  1996]. 


As  mentioned,  the  source  of  these  problems  is  the  non-uniform  coverage  of  the  available  data. 
This  is  a  particular  problem  at  low-latitudes  because  vertical  E  xB  drift  leads  to  the  formation  of 
the  equatorial  anomaly  and  the  development  of  a  complex  morphology  quite  different  from  that 
of  the  mid  latitudes.  However,  the  problems  at  low  latitudes  are  also  exacerbated  by  the  large 
day-to-day  variability  of  this  region.  In  particular,  it  is  problematic  at  best  to  attempt  to  model  a 
feature  such  as  the  equatorial  anomaly  in  terms  of  monthly  medians  when  all  the  characteristics 
of  the  anomaly  (intensity,  location,  extent,  latitudinal  shape)  are  known  to  vary  greatly  day-to- 
day.  For  example,  due  to  the  changing  position  and  extent  of  the  anomaly  a  single  low  latitude 
ground  station  can  day-to-day  effectively  be  sampling  different  geophysical  regions.  Assembling 
data  averages  from  several  such  stations  can  result  in  a  smeared  out  anomaly  that  is  not 
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representative  of  an  actual  anomaly  [Rawer,  1995].  While  for  some  purposes  such  an  average 
ionosphere  is  all  that  is  needed,  for  other  uses  it  is  critical  to  have  a  representative  or  "typical" 
ionosphere. 

One  effect  of  the  issues  described  above  is  that  theoretical  modeling  has  begun  to  appear  as  a 
complementary  approach  to  strict  empirical  modeling.  Rush  et  al.  [1983,  1984]  used  theoretical 
modeling  to  fill  in  the  data  gaps  over  the  mid-latitude  oceans  and  along  with  ground-based 
ionosonde  observations  produced  an  improved  set  of  maps  of  the  F2  peak  critical  frequency 
(f„F2).  Such  maps  are  the  sources  used  by  IRI  for  f0F2  and  hence  NmF2.  Further  improvements  to 
the  f0F2  maps  were  made  using  theory  over  the  oceans,  a  larger  data  base  of  ionosondes,  and  a 
harmonic  analysis  that  was  latitude  dependent  [Fox  and  McNamara,  1988].  At  about  the  same 
time,  efforts  appeared  that  used  first-principles  theoretical  models  to  produce  parameterized 
models  that  would  describe  low  latitude  theoretical  morphology  [Anderson  et  al.,  1987,  1989], 
The  purpose  was  to  produce  a  representative  equatorial  anomaly  rather  than  a  smeared  out 
average  anomaly  [Daniell  et  al.,  1995].  In  this  paper,  we  propose  a  new  option  for  IRI  which 
would  allow  the  results  of  low  latitude  theoretical  modeling  of  electron  density  to  be 
incorporated  directly  into  IRI.  The  purpose  is  to  both  provide  an  improvement  to  IRI  as  well  as 
to  illustrate  one  way  in  which  our  theoretical  understanding  can  be  used  to  fill  in  the  gaps  of  the 
observational  databases  used  by  IRI.  In  the  next  section,  we  present  some  observations  compared 
with  IRI-90  in  order  to  quantitatively  illustrate  some  of  problems  outlined  above.  This  is 
followed  by  a  description  of  the  theoretical  modeling  of  the  low  latitude  ionosphere,  an  example 
of  a  PIM  option  to  IRI,  and  finally  a  discussion  and  conclusion  section. 

Section  2  -  Observations 


In  this  section,  we  originally  presented  a  comparison  between  IRI  and  Jicamarca  observations 
from  a  particular  day.  The  point  of  this  comparison  was  to  simply  illustrate  the  type  of  problem 
that  can  be  encountered  at  low  latitudes  when  comparing  IRI  and  actual  observations.  However, 
given  that  IRI  is  a  model  of  monthly  profiles,  it  should  not  be  expected  to  necessarily  reproduce 
the  profiles  on  a  given  day.  While  the  broad  shape  of  the  afternoon  profiles  and  the  post  sunset 
lifting  of  the  layer  are  not  rare  events  and  occur  regularly,  in  order  to  validate  that  IRI  has  in  fact 
a  low  latitude  problem,  comparisons  need  to  be  made  between  some  sort  of  averaged  data  and 
IRI.  In  Figures  21-23,  which  are  adapted  from  Creamer  [1992],  we  show  the  hourly  averaged 
data  from  several  years  of  Jicamarca  observations  for  moderate  solar  flux  conditions.  The 
December  solstice  refers  to  data  taken  during  November  through  February  and  the  June  solstice 
refers  to  data  from  May  to  August.  The  equinox  period  was  taken  to  be  March  and  April  along 
with  September  and  October.  These  plots  also  include  results  from  IRI-90  and  the  PIM.  In  this 
section,  we  will  focus  on  the  data  (solid  curves)  and  IRI-90  (dashed  curves)  results.  PIM  results 
(dotted  curves)  will  be  discussed  in  the  following  section.  The  IRI  runs  were  made  using  a 
sunspot  number  of  80  and  for  days  of  the  year  1, 100,  and  182.  The  data  are  from  days  where  the 
F10.7  ranged  from  100  to  150. 

We  begin  by  noting  that  the  observed  profiles  (solid  curves)  are  generally  broader  than  IRI-90 
(dashed  curves)  though  in  several  cases  IRI-90  does  a  reasonable  job.  In  Figure  21  (1200  LT)  we 
see  that  for  the  December  solstice,  the  shape  of  IRI's  profile  compares  quite  well  with  the 
averaged  data.  At  equinox,  we  begin  to  see  some  differences  in  both  the  topside  and  bottomside 
shapes  between  IRI  and  data.  For  the  June  solstice,  the  differences  in  the  bottomside  shape  are 
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FIGURE  21.  Density  Profiles  of  PIM,  IRI-90  and  Jicamarca  Data  1200  LT 
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FIGURE  22.  Density  Profiles  of  PIM,  IRI-90  and  Jicamarca  Data  1800  LT 
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FIGURE  23.  Density  Profiles  of  PIM,  IRI-90  and  Jicamarca  Data  2000  LT 
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most  apparent.  In  Figure  22  (1800  LT)  we  see  similar  trends  with  clear  differences  between  IRI- 
90  and  data  in  topside  shapes  during  the  December  solstice  and  the  equinox  and  bottomside 
differences  during  the  June  solstice  and  the  equinox.  Finally,  in  Figure  23  (2000  LT),  the  post 
sunset  lifting  of  the  layer  is  evident  (solid  curves)  during  the  December  solstice  and  the  equinox 
and  is  seen  to  be  missed  by  IRI-90  (dashed  curves).  There  is  a  smaller  amount  of  lifting  of  the 
layer  for  the  June  solstice,  but  it  is  difficult  to  see  given  that  the  data  was  placed  into  40  km  bins 
when  producing  these  average  profiles.  If  hmF2  is  taken  from  the  unbinned  data  and  then 
averaged,  the  June  solstice  results  then  show  a  peak  rise  of  just  16  km  between  1700  LT  and 
2000  LT.  In  contrast,  the  equinox  shows  a  55  km  increase.  Similar  features  are  seen  for 
maximum  solar  flux  conditions  but  not  for  solar  minimum  conditions.  IRI-90  also  consistently 
shows  an  increased  scale  height  above  700  km.  While  the  data  shows  some  evidence  that  this  is 
not  the  case  (see  equinox  at  1200  LT  and  1800  LT),  in  general  these  observations  are  too  low  in 
altitude  to  clearly  determine  whether  this  increased  scale  height  is  appropriate. 

While  the  observations  shown  thus  far  come  from  one  location,  similar  features  are  also  evident 
in  other  longitude  sectors.  Many  of  the  works  referenced  in  the  introduction  describe  studies  that 
included  data  from  other  longitude  sectors.  In  Figure  24,  we  present  profiles  [R.  Tsunoda, 
personal  communication,  1995]  observed  by  the  Kwajalein  incoherent  scatter  radar  located  near 
the  equator  in  the  Pacific  sector.  On  the  left  is  a  profile  from  the  early  morning  before  the 
upward  drift  E  x  B  drift  has  acted  to  broaden  the  peak,  while  on  the  right  we  see  later  in  the  day 
the  characteristic  broad  layer  of  the  low-latitude  daytime  that  results  from  the  upward  E  x  B  drift. 
Further,  it  has  been  clearly  established  that  the  equatorial  anomaly  is  global  in  nature  [Walker, 
1981;  Su  et  al.,  1996]  and  that  the  E  x  B  drift  though  different  in  detail  is  similar  at  all  longitudes 
[Fejer  et  al,  1995].  Given  that  the  discrepancies  described  above  are  related  to  the  effects  of 
vertical  E  x  B  drift  and  the  equatorial  anomaly,  it  is  then  consistent  that  IRI  would  have  these 
low  latitude  problems  at  other  longitudes. 

Section  3  -  Theoretical  Modeling 

In  this  section,  we  expanded  our  discussion  to  clarify  the  relation  between  PIM  and  GTIM  and 
included  comparisons  between  PIM  and  averaged  Jicamarca  data. 

Over  the  last  three  decades,  there  have  been  many  theoretical  studies  of  the  equatorial  F2  region 
and  a  number  of  computer  models  have  been  developed  that  reproduce  the  basic  features  of  this 
region.  Most  recently,  these  models  have  evolved  to  the  point  that  given  realistic  inputs  for  a 
particular  day  they  are  capable  of  reproducing  the  ionosphere  of  that  day  [Preble,  et  al.,  1994]. 
One  such  model  is  the  GTIM  which  calculates  electron  density  profiles  as  a  function  of  location 
and  time  by  solving  ion  continuity  and  momentum  equations.  For  the  F-region,  the  0+  ion 
density  is  determined  by  numerically  solving  the  time-dependent  ion  (0+)  continuity  equation. 
The  needed  velocity  parallel  to  the  magnetic  field  comes  from  combining  the  electron  and  ion 
momentum  equations.  The  perpendicular  velocity  is  taken  as  the  E  x  B  drift  where  the  electric 
field  has  to  be  provided  as  an  input  to  the  model.  A  variety  of  other  geophysical  inputs  are  also 
required  including  the  neutral  constituents,  temperature,  and  winds;  the  electron  and  ion 
temperatures;  the  solar  extreme  ultraviolet  (EUV)  fluxes;  and  the  geomagnetic  field.  For  the 
lower  altitudes,  the  molecular  ions  are  determined  by  numerically  solving  the  steady-state  local 
approximation  of  four  coupled  ion  continuity  equations  [Anderson  et  al.,  1996].  In  Figure  25,  we 
present  vertical  profiles  for  1200  and  1900  LT  from  GTIM  simulations  (dotted  curves)  that  used 
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FIGURE  24.  (a)  Profiles  from  Kwajalein  Radar  Early  Morning 
(b)  Profiles  from  Kwajalein  Radar  Late  in  Day 
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FIGURE  25.  (a)  GTIM  and  Jicamarca  Profiles  1200  LT 
(b)  GTIM  and  Jicamarca  Profiles  1900  LT 
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the  Jicamarca  observed  upward  E  x  B  drift  of  October  1  and  2, 1970.  We  also  include  in  the  plot 
the  electron  density  profiles  (solid  curves)  observed  by  the  Jicamarca  radar  at  those  times.  We 
see  that  when  the  observed  drift  is  used  as  an  input  to  the  model  that  the  agreement  between 
model  and  observation  is  excellent.  Further,  we  see  that  the  model  reproduces  both  the  broad 
profile  shape  of  the  daytime  and  the  lifting  of  the  F2  layer  during  the  post-sunset  period.  The 
additional  curves  (dashed  curves)  are  GTIM  results  when  zero  vertical  E  x  B  drift  was  used  as 
input.  These  simulations  illustrate  the  critical  role  that  the  magnitude  of  the  E  x  B  vertical  drift 
has  in  determining  both  the  shape  of  the  daytime  profile  and  the  post  sunset  rising  of  the  F2  iayer- 

Recently,  a  computationally  fast  PIM  was  developed  at  the  Phillips  Laboratory  Geophysics 
Directorate  as  part  of  a  global,  real-time  ionospheric  specification  model  called  PRISM  which  is 
presently  operational  at  the  Air  Force  50th  Weather  Squadron  (50WS)  in  Colorado  Springs 
[Daniell  et  al.,  1995].  PIM  is  based  on  many  runs  of  several  theoretical  ionospheric  models  and, 
as  such,  is  a  global  ionospheric  model  based  on  theoretical  climatology  rather  than  empirical 
climatology.  For  the  low  latitudes,  the  theoretical  model  used  was  an  early  version  of  GTIM.  It 
was  run  under  three  solar  cycle  conditions  (low,  middle,  and  high  F10.7  cm  flux  values)  and 
three  seasons  (fall/spring  equinox,  summer,  and  winter  solstice  periods)  for  four  longitude 
sectors  (American,  Brazilian,  European/Indian,  and  Pacific).  Electron  density  profiles  every  2° 
latitude  and  every  hour  local  time  were  parameterized  by  9  Empirical  Orthonormal  Functions 
(EOFs).  At  low  and  mid-latitudes  these  EOFs  were  analytically  fit  in  latitude  and  kept  in  tabular 
form  over  the  24-hour  period.  These  tables  and  functions  comprise  the  low-latitude  portion  of 
PIM  and  reproduce  very  accurately  the  electron  density  profiles  generated  by  GTIM. 

Given  the  success  of  theoretical  models  in  reproducing  the  basic  features  of  the  low  latitudes,  we 
would  expect  that  the  PIM  should  likewise  reproduce  those  features.  In  Figures  21-23,  we  have 
included  PIM  results  (dotted  curves)  along  with  the  averaged  Jicamarca  data  (solid  curves)  and 
the  IRI-90  (dashed  curves)  results.  In  Figure  21  for  the  December  solstice,  we  see  that  the  PIM 
produces  a  profile  that  is  broader  than  the  data.  Both  the  PIM  and  IRI-90  underestimate  NmF2 
and  to  some  extent  the  topside  and  bottomside  data  seem  to  split  the  difference  between  PIM  and 
IRI-90.  At  the  equinox,  PIM's  (dotted  curve)  topside  shape  gives  good  agreement  with  data 
(solid  curve)  up  to  around  600  km.  Above  that  height  both  the  PIM  and  IRI  miss  the  sharp  fall 
off  in  the  data.  Again  in  the  bottomside,  the  data  distribution  splits  the  difference  between  the 
models,  though  if  we  normalize  IRI  to  the  F2  peak  of  the  observations,  IRI  (dashed  curve)  has 
the  better  bottomside  shape.  Finally  for  June  solstice,  the  PIM  produces  a  broader  layer  than  IRI- 
90  with  the  data  falling  between  the  two  models.  However,  this  time  it  is  the  PIM  that  has  the 
better  bottomside  shape  with  IRI-90  showing  a  much  thinner  bottomside  than  seen  in  the  data. 
Turning  to  Figure  22  and  1800  LT,  we  again  see  a  broader  PIM  (dotted  curves)  as  compared  to 
IRI-90  (dashed  curves).  For  December,  the  PIM's  topside  agrees  quite  nicely  with  the  data  while 
IRI-90  shows  a  significantly  different  shape.  Moving  to  equinox,  up  to  around  600  km  the  PIM 
gives  reasonable  agreement  with  the  observed  topside  but  as  was  the  case  at  1200  LT  both  the 
PIM  and  IRI-90  miss  the  sharp  fall  off  seen  in  the  data  at  higher  altitudes.  In  the  bottomside  the 
observed  layer  appears  to  be  a  bit  thicker  than  either  model.  The  June  solstice  again  shows  an 
observed  topside  shape  that  falls  between  the  two  models  and  a  observed  bottomside  shape  that 
is  thicker  than  either  model.  At  2000  LT  (Figure  23)  similar  trends  continue.  The  PIM  gives 
good  agreement  with  observations  over  much  of  the  F  region.  In  particular,  the  PIM  captures  the 
lifting  of  the  layer  that  was  missed  by  IRI-90.  For  equinox,  the  PIM  captures  the  topside  shape 
and  the  lifting  of  the  layer.  However,  both  IRI-90  and  PIM  produce  much  thinner  bottomsides 
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than  seen  in  the  averaged  data.  Finally  for  June  solstice,  the  PIM  shows  a  layer  that  has  risen  to  a 
higher  altitude  ( hmF2  ~  440  km)  than  that  seen  in  the  observations  ( h,„F2  ~  410  km),  while  IRI-90 
is  down  at  a  lower  altitude  (hmF2  ~  355  km). 

Having  seen  how  both  IRI  and  the  PIM  compare  to  data  from  Jicamarca  for  solar  moderate 
conditions,  let  us  now  compare  the  two  models  over  a  larger  region  and  for  solar  maximum 
conditions.  In  Figure  26,  we  directly  compare  the  PIM  and  IRI-90  profiles  at  1700  UT  in 
Jicamarca  longitude  meridian.  The  local  time  is  around  1200  and  we  see  that  the  PIM  produces  a 
much  broader  F  layer  around  the  magnetic  equator  than  does  IRI-90.  However,  the  PIM  anomaly 
falls  off  much  faster  with  latitude  than  does  IRI-90  so  that  20°  north  of  the  magnetic  equator 
(around  10°  north  geographic)  IRI-90  is  producing  a  thicker  layer.  This  illustrates  how  in  regions 
of  great  variability  theoretical  models  can  produce  sharper  features  (in  this  case  an  anomaly  that 
is  narrower  in  latitude)  than  may  be  seen  in  models  based  on  averaged  data.  In  Figure  27,  it  is 
h„,F2  that  is  compared  over  the  entire  globe  at  0  UT  and  we  see  that  the  largest  differences 
between  the  PIM  and  IRI-90  are  confined  to  the  post  sunset  period  near  the  equator.  The  topside 
half  thicknesses  are  compared  in  Figure  28  and  here  a  sizable  difference  is  seen  near  the 
magnetic  equator  but  extended  over  180°  in  longitude.  This  also  points  out  how  the  differences 
in  thickness  between  models  are  not  confined  to  just  one  longitude  sector.  In  Figure  26,  a  clear 
difference  in  layer  thickness  was  seen  in  the  American  sector  at  1700  UT.  At  1200  UT,  this 
difference  is  seen  in  the  European/Indian  sector. 

Section  5 


While  Section  4  required  few  changes,  Section  5  was  considerably  expanded  by  including  a 
discussion  of  the  validation  of  the  PIM  as  well  as  the  effect  of  using  IRI-95  versus  IRI-90. 

While  we  have  demonstrated  the  possibility  and  potential  of  the  PIM  option,  the  same  type  of 
validation  is  needed  for  the  PIM  as  is  necessary  for  any  version  of  IRI.  However,  validating  the 
PIM  option  faces  the  same  problems  that  have  hampered  IRI  validations;  that  is,  the  availability 
of  suitable  observations.  To  date,  validation  of  the  low-latitude  portion  of  the  PIM  has  been 
confined  to  the  American  sector,  primarily  to  comparisons  with  data  from  the  Jicamarca 
incoherent  scatter  radar.  Further  validation  for  the  PIM  option  awaits  other  data  sources.  One 
near  term  possibility  is  the  TEC  observations  that  are  available  from  the  dual  frequency 
dispersive  radar  on  the  TOPEX/POSEIDON  satellite  [Bilitza  et  al,  1996].  These  observations! 
which  began  around  the  beginning  of  1993,  would  appear  to  offer  an  exciting  opportunity  for  low 
latitude  testing  of  IRI,  PIM  and  GTIM.  J 

A  key  issue  when  considering  the  validation  of  models  such  as  the  IRI  or  the  PIM  is  that  of  day- 
to-day  variability.  This  is  especially  critical  in  a  region,  such  as  the  low  latitudes,  where  that 
variability  can  be  quite  large.  This  issue  manifests  itself  in  several  ways.  For  example,  the 
process  by  which  one  constructs  an  average  from  data  becomes  critical  especially  if  dealing  with 
any  type  of  variable  structure.  As  mentioned  in  the  introduction,  the  latitude  extent  of  the 
anomaly  can  easily  be  smeared  out  if  one  takes  routine  averages  of  ground  based  measurements. 
Another  example  is  that  of  the  shape  of  the  altitude  profile  changing  markedly  from  day  to  day 
even  when  the  solar  and  magnetic  conditions  are  not  markedly  different.  In  that  case  an  "average 
profile"  may  represent  only  one  of  the  possible  profile  shapes.  In  Figure  21  for  the  December 
solstice,  we  see  an  average  profile  with  a  topside  shape  falling  between  IRI-90  and  the  PIM.  If 
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one  examines  the  scatter  plots  of  the  18  profiles  that  went  into  making  up  this  average  [Creamer, 
1992],  you  see  a  range  of  profiles  that  go  from  very  broad  shapes,  low  peak  densities  (7  x  105), 
and  high  hmF2 s  (>  400  km)  to  thin  topsides  (thinner  than  IRI),  high  peak  densities  (2.5  x  106),  and 
lower  hmF2 s  (~  325  km).  In  this  case,  both  IRI  and  the  PIM  fall  within  this  range  of  profile 
shapes.  We  can  speculate  that  the  PIM  profile  shape  is  representative  of  those  days  when  the 
electric  field  was  of  a  certain  strength,  while  IRI,  based  on  solar  minimum  data,  possibly 
represents  those  profiles  from  days  with  weaker  electric  fields.  But  all  the  cases  fall  into  the 
same  category  of  magnetically  quiet,  moderate  solar  conditions.  Of  course  with  such  a  range  of 
shapes,  it  is  difficult  to  determine  what  is  "typical"  and  there  is  little  reason  to  expect  that  a 
theoretical  model  driven  by  an  average  E  x  B  drift  based  on  one  limited  set  of  drift  data  will 
necessarily  produce  the  "average"  seen  in  a  limited  set  of  electron  density  profile  data  from  a 
different  period.  In  such  regions  with  large  variability,  what  is  needed  for  IRI  is  an  expanded 
capability  that  quantifies  the  variations  about  the  "average  profile."  As  for  the  PIM,  it  is  planned 
that  the  next  version  will  be  parameterized  by  some  sort  of  electric  field  strength  parameter.  In 
order  to  validate  such  a  version  of  the  PIM  and  the  above  speculation,  it  would  be  useful  to 
perform  a  GTIM  model  study  where  a  database  of  simultaneous  vertical  drift  and  electron  density 
profile  measurements  are  available.  In  that  case,  one  could  verify  whether  the  observed  day-to- 
day  variations  in  profile  shape  can  be  modeled  as  directly  coming  from  the  observed  day-to-day 
variations  in  vertical  drift. 

The  scale  height  differences  above  700  km  between  the  IRI-90  and  the  PIM  profiles  are  an 
obvious  issue  that  needs  to  be  addressed  as  we  consider  the  issue  of  model  validation.  As 
mentioned  in  Section  2,  the  increased  scale  height  shown  in  IRI-90  suggests  a  transition  from  an 
O  dominated  F  region  to  an  H+  dominated  topside.  This  suggestion  comes  from  both  theory  and 
observations  that  have  shown  that  this  type  scale  height  change  occurs  when  this  transition  of 
species  occurs  [Bilitza,  1990],  Further,  the  transition  heights  observed  at  solar  minimum  are 
generally  consistent  with  the  transition  heights  evident  in  the  IRI-90  profiles.  That  this  is  so  is  of 
course  reasonable  given  that  IRI's  topside  is  based  on  solar  minimum  topside  sounder  data. 
However,  observations  have  shown  that  the  transition  height  moves  up  in  altitude  at  higher  levels 
of  solar  activity.  This  has  been  verified  recently  in  a  study  of  DMSP  satellite  data  [West  et  al., 
1997]  that  showed  during  low  solar  activity  the  transition  height  is  well  below  800  km  and  well 
above  800  km  during  higher  solar  activity.  Further  for  the  times  of  the  DMSP  data,  it  was  seen 
that  the  W  density  at  800  km  never  gets  much  above  1  or  2  x  104.  In  disagreement  with  this  are 
the  IRI-90  results  that  show  for  solar  moderate  conditions  at  21.8  LT  (the  time  of  DMSP  data)  a 
scale  height  change  at  around  700  km  and  H*  densities  at  800  km  of  around  9  x  104.  An 
observational  approach  to  this  problem  should  come  from  topside  sounder  data.  However,  only  a 
small  part  of  the  Alouette-ISIS  topside  sounder  database  has  been  available  for  use  [Bilitza, 
1994].  It  is  expected  that  if  and  when  more  of  this  database  becomes  available  significant 
improvements  can  be  made  in  IRI's  topside  morphology.  On  the  theoretical  side,  we  find  that  the 
PIM  does  not  include  H*  and  so  obviously  cannot  describe  the  transition  between  species. 
Fortunately,  a  new  version  of  GTIM  is  under  development  that  includes  I T  at  low-and  mid¬ 
latitudes.  It  is  anticipated  that  the  next  version  of  the  PIM  will  be  produced  using  this  improved 
GTIM  and  thus  should  include  the  transition  from  0+  to  IT  that  occurs  in  the  topside  The  two 
issues  discussed  above  (daily  profile  shape  changes  and  OW  transitions)  are  examples  of  the 
problems  that  arise  in  part  from  the  lack  of  suitable  databases.  They  are  also  examples  of  how 
theoretical  models  can  be  used  to  complement  pure  empirical  modeling.  That  is,  validated  theoiy 
can  serve  to  fill  in  the  data  gaps  and  provide  a  spatial  and  temporal  framework  when  databases 
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are  too  sparse  to  do  so.  We  expect  that  a  next  generation  PIM  that  includes  both  tF  and  an 
electric  field  parameterization  can  play  such  a  role  and  serve  to  complement  strict  empirical 
modeling. 

Finally,  this  work  was  done  using  IRI-90;  however,  during  the  writing  of  this  paper  IRI-95 
became  available  and  could  be  run  interactively  via  a  World  Wide  Web  site.  We  have  rerun  the 
cases  shown  in  Figures  21-23  using  IRI-95  and  find  essentially  no  differences  in  the  topside  and 
F?  peak  parameters.  However  for  day  1  and  100,  the  bottomsides  at  12  LT  coming  from  ERI-95 
are  somewhat  thinner  than  the  IRI-90  results.  Given  these  results,  it  is  evident  that  essentially  all 
our  conclusions  based  on  IRI-90  remain  true  for  IRI-95. 

2.5.  Ionospheric  Modeling  Using  Systems  of  Plasma  Kinetic  Equations 

In  certain  regions  of  the  ionosphere,  plasma  kinetic  equations  can  be  useful  for  studying 
problems  where  wave-particle  interactions  are  important  sources  of  particle  heating  and 
transport.  However,  a  general  system  of  kinetic  equations  coupled  with  Maxwell's  equations 
forms  a  formidable  set  of  equations.  One  approach  for  solving  this  set  is  to  assume  that  the 
fluctuating  E  and  B  fields  are  known  and  can  be  used  in  a  physical  model  that  describes  the 
wave-particle  interactions.  This  model  is  then  used  in  the  set  of  equations  to  solve  for  the  one- 
particle  distribution  function  and  the  quasistatic  electric  field.  However,  difficulties  still  remain 
in  self-consistently  solving  for  the  quasistatic  (ambipolar  field)  electric  field  and  the  distribution 
function.  In  particular,  it  is  not  practical  to  attempt  a  direct  numerical  solution  of  Poisson's 
equation.  The  approach  being  explored  here  is  to  use  a  multi-energy-group  (MEG)  method  to 
derive  moment  equations  and  to  solve  these  equations  for  the  moments  and  the  ambipolar  field. 
This  electric  field  can  then  be  used  in  solving  a  kinetic  equation  for  the  one-particle  distribution 
function. 

A  multi-energy-group  method  of  analyzing  a  system  of  plasma  kinetic  equations  involves 
defining  the  one-particle  distribution  function  for  each  type  of  particle  as  a  sum  of  functions,  one 
for  each  major  energy  group,  and  then  constructing  a  set  of  moment  equations.  The  set  of 
equations  that  then  need  to  be  solved  depend  on  the  number  of  energy  groups  chosen  and  the 
method  used  to  close  the  set  of  moment  equations.  We  have  been  working  on  two  versions  of 
this  method.  One  produces  a  set  of  8  differential  equations  while  the  second  involves  11 
equations.  These  equations  describe  how  8  or  1 1  variables  vary  in  one  spatial  dimension.  A 
computer  program  was  written  that  performs  the  symbolic  manipulation  needed  to  transform 
these  sets  of  equations  into  a  triangular  form  so  they  can  be  easily  solved  using  standard 
techniques.  This  program  also  calculates  the  necessary  derivatives  and  generates  the  FORTRAN 
code  that  is  needed  for  the  subroutines  that  are  used  by  a  standard  differential  equation  solver.  A 
driver  program  that  uses  the  equation  solver  has  been  written  and  the  first  preliminary  solutions 
for  these  sets  of  equations  have  been  obtained. 

Work  began  on  developing  a  numerical  solution  to  the  kinetic  equation.  The  first  attempts  have 
involved  solving  a  differential  equation  in  two  velocity  variables  and  one  spatial  variable.  A  split 
operator  method  is  used  with  one  operator  describing  the  effects  of  gravity,  the  ambipolar 
electric  field,  and  the  magnetic  field  (mirroring  force)  while  the  other  describes  the  effects  of  the 
wave-particle  interactions.  The  most  difficult  step  in  writing  the  algorithm  has  been  the  proper 
handling  of  the  coordinate  grids  as  these  have  to  be  modified  every  half  step  as  the  equation  is 
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integrated.  A  preliminary  code  was  implemented  and  tested  on  several  simple  examples  of  the 
differential  equation.  The  preliminary  code  was  then  rewritten  several  times  to  improve  the  way 
the  coordinate  grids  were  varied  at  each  half  step.  In  addition,  numerical  integration  of  results  at 
the  full  step  was  used  to  see  if  the  various  moment  equations  were  being  satisfied.  In  several 
preliminary  runs  of  the  latest  code,  it  was  found  that  the  moment  equations  were  not  being 
satisfied.  Work  continues  to  identify  the  reason  for  this. 


3.  OBSERVATIONAL  DATABASES 

3.1.  Solar  Cycle  Dependence  of  Absolute  Ionospheric  Total  Electron  Content 

The  aim  of  this  research  is  to  aid  in  the  development  of  an  improved  single  frequency  GPS  user 
algorithm  by  defining  an  improved  fit  between  TEC  and  long-term  solar  activity  and  to 
determine  if  that  fit  changes  for  different  seasons  and  geographic  locations.  The  ionosphere  is 
the  most  substantial  source  of  error  in  GPS  positioning  and  navigation.  Since  the  ionosphere  is  a 
dispersive  medium,  the  dual-frequency  GPS  user  can  automatically  correct  for  ionospheric 
effects.  For  the  single  frequency  GPS  user,  the  Ionospheric  Correction  Algorithm  OCA)  was 
developed  to  correct  for  approximately  50  percent  rms  of  the  ionospheric  error.  The  ICA  was 
developed  in  1975  using  the  Bent  Model,  a  state-of-the-art  empirical  ionospheric  model  that  was 
based  on  data  recorded  during  an  average  solar  cycle.  Solar  cycle  maximum  values  since  that 
period  have  been  much  higher.  As  a  result,  the  ICA  does  not  adequately  represent  periods  of 
extremely  high  solar  maximum  activity.  In  addition,  the  ICA  is  deficient  in  describing  the 
significant  ionospheric  changes  that  are  characteristic  of  the  equatorial  and  low-latitude  regions. 
The  study  is  accomplished  by  correlating  values  of  the  10.7  cm  solar  radio  flux  (F107)  with  a 
worldwide  data  base  of  foF2  measurements  recorded  from  1958  through  1992,  a  period  covering 
almost  four  solar  cycles.  This  study  is  further  augmented  with  TEC  from  Faraday  rotation 
measurements  and  TEC  from  dual-frequency  GPS  data  recorded  at  Goldstone  CA.  A  special 
emphasis  will  be  placed  on  resolving  whether  there  is  a  real  saturation  effect  in  TEC  at  extreme 
levels  of  10.7  cm  solar  flux  as  evidenced  in  work  by  Balan  et  al.  [1993, 1994]. 

The  results  of  this  study  illustrate  that  day-to-day  variations  in  Fi07  do  not  correlate  well  with 
day-to-day  changes  in  TEC.  In  particular,  one  day  spikes  in  F10.7  do  not  result  in  spikes  in  TEC. 
This  is  likely  because  spikes  in  solar  radio  flux  are  not  well  correlated  with  spikes  in  solar  EUV 
radiation,  the  primary  source  of  ionization  in  the  F2  region  of  the  ionosphere.  This  study  further 
illustrates  that  a  much  better  correlation  between  solar  activity  and  observed  TEC  is  obtained  by 
using  a  27  day  running  mean  of  F10.7  rather  that  the  daily  value  or  shorter  period  (e.g.,  5  day) 
running  means.  This  smoothes  out  the  individual  spikes  and  results  in  better  correlation  between 
solar  activity  and  observed  TEC.  This  result  is  important  in  the  development  of  an  improved 
GPS  single  frequency  correction  algorithm. 

As  part  of  this  study,  we  have  continued  to  access  and  analyze  dual  frequency  GPS 
measurements  from  the  GPS  International  GPS  Geodynamics  Service  for  (IGS)  maintained  by 
the  Jet  Propulsion  Laboratory.  As  of  the  end  of  this  year,  we  have  completed  data  processing  for 
the  data  recorded  through  August  1996  at  four  different  geographic  locations.  In  addition,  the 
GPS  measurements  made  at  Goldstone  CA  from  December  90  through  June  96  have  been 
utilized  to  estimate  the  potential  errors  made  with  the  current  GPS  single  frequency  algorithm. 
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Unfortunately,  there  is  a  significant  gap  in  the  data  from  September  91  through  July  92.  This  gap 
limits  our  attempts  to  quantify  the  potential  errors  of  the  current  single  frequency  algorithm 
during  solar  maximum.  Figure  29  illustrates  the  results  of  the  Goldstone  analysis.  In  the  top 
three  figures,  daily  values  of  10.7  cm  solar  flux  are  plotted  versus  dual-frequency  GPS 
measurements  made  at  1400  hours  local  time  for  each  of  three  seasons,  winter  (November, 
December,  January  and  February),  equinox  (March,  April,  September  and  October)  and  summer 
(May,  June,  July  and  August).  The  bottom  three  figures  illustrate  the  corresponding  predictions 
of  ionospheric  behavior  for  the  Goldstone  site  from  the  current  single-frequency  ionospheric 
correction  algorithm.  In  comparing  the  measured  data  with  the  algorithm  predictions,  one  notes  a 
similar  linear  trend  up  to  an  F10.7  value  of  approximately  210.  At  that  point,  the  ionospheric 
correction  algorithm  underestimates  the  measured  data.  Much  of  this  apparent  saturation  effect 
is  caused  by  the  strong  one  day-peaks  in  F10.7. 

3.2.  Absolute  Real-Time  Ionospheric  Measurements  from  GPS  Satellites  in  the  Presence 
of  Anti-Spoofing 

GPS  satellites  transmit  coherent  radio  signals  at  two  frequencies,  LI  at  1.575  GHz  and  L2  at 
1.227  GHz.  These  carriers  are  phase-modulated  with  two  types  of  code,  coarse  acquisition  (C/A) 
on  LI  and  precise  code  (P-code)  on  LI  and  L2.  The  P  codes  are  used  to  determine  the  group 
delay  measurements.  Absolute  measurements  of  the  ionospheric  range  delay  can  be  obtained  by 
differencing  the  group  delay  from  both  frequencies.  On  31  January  1994,  anti-spoofing  (a/s)  was 
turned  on  for  most  block  II  GPS  satellites.  When  a/s  is  turned  on  the  precision  code  (P-code)  is 
encrypted  and  replaced  by  an  unknown  Y-code.  Only  military  receivers  equipped  with  an 
authorized  key  can  use  the  Y-code  directly.  In  this  study,  we  have  accessed  dual-frequency  GPS 
data  recorded  with  GPS  Rogue  and  Turborogue  receivers  designed  by  Allen  Osborne  Associates, 
Inc..  When  a/s  is  activated,  these  receivers  continue  to  provide  group  delay  and  carrier  phase 
measurements  in  a  codeless  mode  by  cross  correlating  the  LI  and  L2  signal.  Because  of  the 
degraded  signal-to-noise  ratio  obtained  by  the  receiver  when  the  satellites  are  in  a/s  mode,  there 
is  a  significantly  higher  level  of  system  noise,  particularly  at  low-elevation  angles.  This 
increased  noise  can  seriously  affect  many  of  the  real-time  applications  for  GPS  that  require  high 
precision  and  accuracy,  such  as  the  precision  aircraft  approaches  in  the  FAA’s  Wide  Area 
Augmentation  System  (WAAS). 

The  purpose  of  this  research  was  to  illustrate  the  limitations  to  which  ionospheric  range  delay  can 
be  measured  in  real-time  under  both  a/s  on  and  off  conditions.  This  was  carried  out  by  testing 
various  filter  techniques  on  real-time  dual  frequency  GPS  data  recorded  at  several  sites  under  a/s 
on  and  off  conditions.  A  summary  of  the  results  is  given  in  the  abstract  of  Gendron  et  al. 
included  in  Section  5,  Presentations  and  Proceedings. 


3.3.  Total  Electron  Content  of  the  Earth's  Protonosphere 

Continuous  measurements  of  the  protonospheric  contribution  to  TEC  from  the  ATS-6  satellite 
were  made  from  a  number  of  stations  covering  middle  and  near-equatorial  latitude  regions  from 
July  1974  through  December  1978,  a  period  of  low  to  moderate  solar  activity.  These  ATS-6 
measurements  are  unique  in  that  they  afforded  the  only  opportunity  to  continuously  monitor  TEC 
using  both  the  Faraday  rotation  and  the  differential  group  delay  technique  of  radio  waves  from 
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FIGURE  29.  (top)  Daily  Values  of  10.7  cm  Solar  Flux  vs  Dual-Frequency  GPS  Measurements  1400  LT 

(bottom)  Predictions  of  Ionospheric  Behavior  from  Current  Single-Frequency  Ionospheric  Correction  Algorithm 


the  same  satellite.  This  combination  of  methods  allowed  the  measurement  of  the  contribution  of 
TEC  from  the  earth’s  protonosphere,  that  region  extending  from  approximately  2,000  km  altitude 
to  the  height  of  the  geostationary  ATS-6  satellite  at  35,800  km.  These  measurements  form  the 
major  source  of  our  knowledge  of  the  geographic  and  temporal  variations  in  the  contribution  to 
TEC  from  the  earth’s  protonosphere.  In  this  quarter,  a  study  was  initiated  to  compare  the 
available  protonospheric  electron  content  observations  from  a  few  representative  stations  with  a 
model  of  the  earth’s  protonosphere.  The  model  will  then  be  used  to  infer  the  protonospheric 
electron  content  contribution  to  TEC  measured  by  the  currently  used  technique  of  differential 
group  delay  from  the  dual -frequency  transmissions  from  the  GPS  satellites. 

Initial  steps  of  this  research  included  a  review  of  some  of  the  more  important  results  of  the  ATS- 
6  experiment  and  some  additional  data  processing  and  analysis  of  the  measurements  recorded  at 
Palehua,  Hawaii.  Monthly  averages  of  TEC  and  protonospheric  electron  content  (Np)  from 
September  1977  through  December  1978  are  shown  in  Figure  30.  In  winter,  the  protonospheric 
content  NP  is  essentially  constant  throughout  the  day  and  night.  In  summer,  there  are  daytime 
maxima  in  the  early  afternoon  (1400  LT).  In  Figure  31,  the  ratio  of  NP  to  TEC  changes  markedly 
with  local  time  with  a  minimum  of  about  10  percent  in  the  early  afternoon  and  a  distinct 
maximum  (between  50  percent  and  90  percent)  near  0500  LT  in  all  months. 

3.4.  Instrumental  Bias  in  GPS  Measurements 

Line-of-sight  ionospheric  measurements  derived  from  differencing  dual-frequency  GPS  range 
data  are  corrupted  by  instrumental  biases  in  both  the  receiver  and  GPS  satellite  transmitters.  The 
instrumental  bias  is  the  difference  of  the  two  dispersive  delays  introduced  by  the  analog  hardware 
in  the  LI  and  L2  signal  paths.  These  instrumental  biases  have  continued  to  be  a  major  concern 
for  ionospheric  research  using  GPS.  Since  it  is  highly  desirable  to  be  able  to  make  measurements 
of  TEC  from  GPS  to  an  absolute  accuracy  of  1  TEC  unit,  (1016  el/m2),  it  is  important  to  be  able  to 
estimate  the  satellite  and  receiver  biases  with  a  very  high  degree  of  accuracy. 

Several  research  groups  have  reported  satellite  and  receiver  bias  values  derived  from  various 
estimation  strategies.  Wilson  and  Mannucci  [1994]  have  evolved  a  technique  that  takes 
advantage  of  a  global  network  of  more  than  50  dual-frequency  stations.  Sardon  et  al.  [1994] 
have  developed  a  technique  that  uses  a  local  network  of  stations  and  generates  results  consistent 
with  Wilson  and  Mannucci  [1994].  These  studies  have  illustrated  that  the  satellite  biases  have 
magnitudes  of  +/-5  nanoseconds  (ns)  of  differential  delay  (1  ns  of  differential  delay  at  L  band  = 
2.85  TEC  units)  and  remain  constant  over  a  fairly  long  period  of  time.  These  studies  have  also 
generated  receiver  biases  for  a  number  of  stations.  Receiver  biases  have  been  shown  to  be  as 
large  as  20  ns  (nearly  60  TEC  units). 

In  our  research  activities  using  GPS,  we  have  frequently  applied  the  satellite  and  receiver  biases 
provided  by  the  Jet  Propulsion  Laboratory  (JPL)  and  described  in  Wilson  and  Mannucci  [1994]. 
Our  experience  has  shown  that  the  satellite  biases  are  consistent  and  dependable.  The  receiver 
biases,  however,  are  not  reliable.  This  is  due  to  the  fact  the  any  change  in  the  receiver,  pre-amp, 
cabling  and  antenna  position  or  configuration  has  a  significant  effect  on  the  overall  receiver  bias 
value.  Unfortunately,  JPL  calculates  receiver  biases  only  on  an  occasional  basis  and  they  are  not 
reported  in  the  data  base.  As  a  result,  the  ionospheric  researcher  is  often  faced  with  guessing  at 
the  receiver  bias. 
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FIGURE  30.  Monthly  Averages  of  TEC  and  Protonospheric  Electron  Content  Sept  77  -  Dec  78 


We  have  attempted  to  validate  a  new  technique  that  provides  estimates  of  combined  satellite  and 
receiver  biases  at  individual  locations.  The  Self-Calibration  of  Pseudorange  Errors  (SCORE) 
[Bishop  et  al.,  1996]  concept  uses  a  self-consistency  constraint  on  the  receiver’s  own 
measurements  of  ionospheric  delay  to  derive  the  sum  of  the  receiver  system  and  pseudorange 
error  for  each  satellite.  The  self-consistency  concept  of  SCORE  is  defined  as  a  constraint  that  the 
ionosphere  changes  only  moderately  over  space  and  time.  Thus,  two  different  satellites  should 
exhibit  a  common  TEC  induced  range  delay  if  their  ray  paths  intersect  the  ionosphere  at  “near 
space-time.”  Other  differences  should  then  be  attributed  to  the  receiver  and  satellite  hardware 
biases.  SCORE  quantifies  “near  space  time”  in  a  weighted  least  squares  approach,  where 
divergence  in  TEC  measurements  are  down-weighted  according  to  their  space-time  separation. 
The  advantages  that  this  method  has  over  those  previously  described  is  that  it  requires  data  from 
only  one  site  and  it  is  faster  computationally.  If  sufficiently  correct,  this  method  may  enable  the 
ionospheric  researcher  to  make  single  site  ionospheric  measurements  at  high  levels  of  accuracy. 

The  approach  of  this  initial  validation  study  was  to  apply  the  SCORE  algorithm  to  dual- 
frequency  GPS  measurements  collected  at  various  mid-and  low-latitude  locations  for  a  series  of 
geomagnetically  quiet  days.  The  results  were  then  compared  for  validity  and  consistency.  The 

data  used  in  this  study  was  recorded  from  16  to  20  September  1995  at  the  following  six 
locations:  c 


Westford,  MA 

42.42°N 

71.49°W 

Bermuda 

32.30°N 

64.70°W 

Richmond,  FL 

25.46°N 

80.38°W 

Bogota,  Columbia 

4.6 1°N 

74.08°W 

Arequipa,  Peru 

16.36°S 

71.49°W 

Santiago,  Chile 

32.97°S 

70.67°W 

These  sites  were  selected  because  they  span  a  wide  latitudinal  range  over  the  American  Sector. 
In  particular,  the  near  equatorial  sites  should  exhibit  extreme  gradients  in  the  ionosphere  and  will 
represent  a  worst  case  scenario  for  SCORE.  Figure  32  summarizes  the  results  of  the  study.  In 
this  figure,  the  combined  satellite  and  receiver  biases  derived  with  SCORE  over  the  5-day  period 
together  with  their  standard  deviations  are  plotted  on  the  y-axis  for  each  of  the  six  recording 
stations.  Note  that  the  biases  have  been  normalized  to  zero  to  enable  a  visual  comparison 
between  the  individual  stations.  In  this  figure,  the  results  exhibit  a  similar  pattern  at  the  three 
mid-latitude  stations  at  Westford,  Bermuda  and  Richmond.  The  Westford  station  illustrates  the 
most  consistent  results  with  very  small  deviations  from  the  mean.  The  Bermuda  and  Richmond 
plots  show  larger  standard  deviations,  although  they  are  still  less  than  3  TEC  units 
(approximately  1  meter).  Results  for  the  near  equatorial  location  at  Bogota  exhibit  very  large 
error  bars  and  the  pattern  seems  different  than  the  one  seen  at  the  mid-latitude  stations.  Results 
at  Arequipa  are  slightly  better  than  those  at  Bogota  but  the  error  bars  are  still  as  large  as  9  or  10 
TEC  units  (nearly  3  meters).  At  Santiago,  the  standard  deviations  are  smaller  but  there  are  large 
discrepancies  in  the  magnitude  of  some  of  the  biases  with  those  illustrated  in  the  northern  mid¬ 
latitude  sector. 

Figures  33  and  34  show  data  analysis  results  when  the  SCORE  biases  that  are  represented  in 
Figure  32  are  utilized.  In  general,  if  the  satellite  and  receiver  bias  calculations  are  correct,  a 
reasonable  diurnal  plot  of  TEC  will  be  obtained  with  high  elevation  measurements.  In  Figure  33, 
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FIGURE  32.  TEC  Values  Using  SCORE  Biases  and  Their  Standard  Deviations 


west  (  42.42,  -71.49)  9/16/95  west 


331  TV3UH3A  !N3TV/\in03  3H  TVOLUOA  1N31WUT03 


46 


BOGT  (  4.61,  -74.08)  9/16/95  BOGT  (  4.61,  -74.08)  9/17/95  BOGT  (  4.61,  -74.08)  9/18/95 


FIGURE  34.  Ionospheric  Effects  At  Bogota 


high  elevation  data  recorded  at  Westford  over  the  5-day  period  is  plotted  using  the  SCORE 
biases.  Note  that  a  reasonable  diurnal  curve  is  achieved  on  all  five  days.  Similar  figures  were 
obtained  with  the  Bermuda  and  Richmond  data.  Figure  34  illustrates  the  data  analysis  results  at 
Bogota.  This  station  is  situated  near  the  equatorial  anomaly  region,  a  region  characterized  by  the 
highest  values  of  TEC  and  peak  electron  density  in  the  worldwide  ionosphere.  In  addition,  the 
day-to-day  variability  of  the  near-equatorial  ionosphere  can  be  very  large,  due  to  changes  in  the 
strength  of  the  E  x  B  drift  at  the  magnetic  equator  and  from  the  differences  in  the  day-to-day 
strength  of  the  meridional  and  zonal  wind  component.  The  results  shown  in  Figure  34  reflect 
these  ionospheric  effects  at  Bogota.  Although  reasonable  diurnal  curves  for  Bogota  are  achieved 
with  the  SCORE  calculations,  the  negative  TEC  values  shown  on  9/18/95  indicate  that  there  are 
errors  in  the  overall  magnitude  of  the  biases. 

The  results  illustrated  in  these  figures  indicate  that  SCORE  provides  excellent  estimation  of  the 
combined  satellite  and  receiver  biases  at  mid-latitude  locations.  Calculations  at  low-latitude 
locations,  however,  are  still  questionable.  Bishop  (1996)  suggests  that  some  of  the  problems 
with  the  SCORE  calculations  at  the  Bogota  and  Arequipa  sites  may  be  related  to  the  large 
ionospheric  gradients  in  the  equatorial  region  but  the  major  problems  are  more  likely  due  to  the 
fact  the  many  of  the  measurements  at  these  sites  were  not  continuous.  He  also  suggested  that  the 
large  day-to-day  deviations  in  the  combined  biases  could  be  related  to  different  satellite  coverage 
and  multi-path  errors  at  the  individual  stations.  In  order  to  affirm  the  performance  of  SCORE  in 
the  low-latitude  region,  it  is  necessary  to  continue  this  study  with  additional  measurements  and  to 
compare  the  results  with  calculations  from  other  bias  estimation  algorithms. 

3.5.  Statistics  of  the  Spatial  and  Temporal  Variations  in  Ionospheric  Range  Delay 

Work  has  been  initiated  to  observe  the  spatial  and  temporal  variations  in  ionospheric  range  delay. 
This  study  will  be  based  on  dual-frequency  GPS  measurements  recorded  from  1993  through  early 
1997  by  the  IGS  network,  managed  by  the  JPL.  GPS  measurements  provide  the  unique 
opportunity  to  observe  simultaneous  ionospheric  estimates  in  up  to  eight  different  slant  paths 
from  a  single  receiver.  These  measurements  can  imply  the  spatial  variations  viewed  from  one 
location.  In  addition,  the  temporal  variation  in  range  delay  at  a  single  location  can  be  made  by 

determining  the  rate  of  change  in  the  ionospheric  as  measured  with  the  highest  elevation  satellite 
at  particular  times  of  day. 

This  analysis  is  intended  to  provide  a  characterization  of  the  ionospheric  variations  that  will  be 
encountered  in  the  Federal  Aviation  Administration’s  (FAA)  operational  WAAS.  The  WAAS 
system  is  expected  to  experience  the  most  problems  during  periods  of  ionospheric  scintillation 
and  during  periods  of  major  geomagnetic  storms  activity  when  the  spatial  and  temporal  gradients 
in  the  ionosphere  can  be  higher  than  usual.  This  study  is  to  be  performed  to  quantify  the 
magnitude  and  frequency  of  these  large  spatial  and  temporal  gradients. 

To  date,  we  have  obtained  and  processed  GPS  data  for  March,  June  and  September  of  1993  and 
for  Januaiy  and  March  of  1994  at  the  following  sites: 

Westford,  Massachusetts 
Goldstone,  California 
Albert  Head,  British  Columbia 
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Fairbanks,  Alaska 
Kokee,  Hawaii. 


Data  analysis  and  statistical  routines  are  continuing  into  the  next  year. 


3.6.  Airglow  Data 

Programs  for  the  SUN  workstation  were  modified  to  accept  the  new  version  of  the  AURIC  Code 
(version  46)  that  was  received.  Considerable  time  was  spent  in  examining  the  AURIC  46  code  in 
anticipation  of  using  it  to  analyze  STS-4,  ST4-39,  and  Polar  Bear  airglow  data.  We  have  also 
examined  programs  that  were  written  for  analysis  of  Polar  Bear  data  by  Bob  Rosenberg  and 
students  at  Tel  Aviv  University,  such  as  converting  the  raw  data  to  images,  filtering  the  data,  roll 
adjustment,  superimposing  the  landmass  on  the  images,  and  showing  the  data  in  the  corrected 
geomagnetic  coordinate  system. 

We  started  work  on  analyzing  data  from  the  STS-4  flight.  Programs  for  plotting  tangent  altitude 
data  and  spectral  data  obtained  on  the  STS-4  flight  were  examined  and  modified  as  needed.  We 
were  then  in  a  position  to  plot  all  the  STS-4  data.  Programs  for  obtaining  tangent  altitude  models 
using  AURIC  for  the  same  parameter  as  the  data  such  as  solar  zenith  angle,  latitude,  longitude, 
Fio.7,  etc.  were  examined.  Several  of  these  models  were  then  plotted  and  compared  with  the  data 
plots.  In  some  cases,  the  data  and  models  agree  very  well,  in  others  only  fair.  Only  a  few  such 
comparisons  have  been  made.  It  is  planned  to  make  the  comparisons  for  all  the  STS-4  tangent 
altitude  data. 

In  order  to  analyze  the  data  obtained  on  STS-39,  we  also  needed  the  ephemeris  data  for  this 
flight.  With  the  cooperation  of  Paul  Pruneau  of  Boston  College,  we  obtained  the  following 
parameters:  time  of  flight  (month,  day,  year,  hour,  minute,  second),  latitude,  longitude,  solar 
zenith  angle,  pitch,  roll,  and  yaw.  We  have  written  programs  to  show  the  pertinent  parameters  on 
the  graphs  depicting  the  intensities  of  the  radiations  measured.  Using  these  parameters,  models 
using  AURIC  were  prepared.  In  doing  this,  averaging  of  spectra  up  to  45,  was  accomplished.  A 
comparison  of  the  models  with  the  data  examined  so  far,  in  the  wavelength  range  of  1100  to 
1700  Angstroms,  showed  good  results  with  the  one  exception  being  the  comparison  at  Lyman 
alpha. 


3.7.  Signature  of  Equatorial  Scintillation  Onset  in  Near-Sunset  DMSP  Satellite  Data 

During  this  year,  we  worked  on  a  study  of  ion  density  data  from  the  DMSP  SSIES  sensor. 
Previous  work  has  shown  that  the  overall  shape  of  latitude  profiles  of  total  ion  density  can  be 
related  to  the  strength  of  the  ionospheric  eastward  electric  field,  and  to  the  strength  and  direction 
of  the  meridional  neutral  wind.  Since  the  occurrence  of  equatorial  spread-F  (a  source  of 
equatorial  radio  signal  scintillation)  is  dependent  on  both  of  these  physical  parameters  it  should 
be  possible  to  correlate  them  with  the  occurrence  of  spread-F.  A  comparison  of  ground-observed 
scintillation  data  from  a  receiver  on  Ascension  Island  (346°  E)  with  corresponding  DMSP  passes 
from  1989  shows  that  days  when  spread-F  did  occur  match  up  well  with  days  when  DMSP 
latitude  profiles  of  ion  density  showed  signatures  of  strong  ionospheric  electric  fields  and  weak 
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neutral  winds.  Further  work  in  this  study  will  compare  scintillation  occurrence  with  DMSP 
satellite  parameters  for  other  longitudes  and  for  different  times  in  the  solar  activity  cycle. 
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5.  PRESENTATIONS  AND  PROCEEDINGS 


We  were  involved  in  15  presentations  at  various  scientific  meetings  of  which  7  appeared  in 
published  proceedings.  There  are  also  4  presentations  that  have  been  submitted  to  meetings  that 
are  to  be  held  in  the  summer  of  1997. 


•  B.  Basu,  D.T.  Decker,  and  J.R.  Jasperse,  “Incident  Proton  Spectra:  Ionospheric  Effects  of 
High  Energy  Power  Law  Tails”,  presented  at  the  Ionospheric  Effects  Symposium  held  in 

345X19%a'  ^  ^  May'  1996  11  a,S°  appeared  in  ]996  ionospheric  Effects  Symposium,  340- 

Abstract 

Studies  of  ion  populations  in  the  central  plasma  sheet  (CPS)  [Christ on  et  al,  1988,  1989, 
1991]  have  found  that  at  high  energies  (E  >  characteristic  or  peak  energy)  there  is  a  nonthermal 
power  law  tail  which  can  be  fitted  using  a  kappa  distribution.  While  similar  studies  for 
ionospheric  altitudes  are  fewer  in  number,  Lyons  and  Evans  [1984]  combined  data  from  various 
instruments  to  show  that  proton  distributions  at  ionospheric  altitudes  have  high-energy  tails  and 
are  the  same  as  those  seen  for  earthward-streaming  protons  in  the  outer  boundary  of  the  plasma 
sheet.  More  recently,  the  Particle  Environment  Monitor  (PEM)  experiment  on  board  the  Upper 
Atmosphere  Research  Satellite  (UARS)  [Sharber  et  al,  1993]  has  observed  ionospheric  ion 
spectra  with  a  high-energy  power  law  tail  similar  to  that  described  by  Christon  et  al.  [1991]  as 
typical  of  CPS  ion  populations.  However,  previous  theoretical  calculations  of  proton-hydrogen 
atom  transport  have  generally  assumed  a  Maxwellian  distribution  for  the  incident  proton  spectra 
[Strickland  et  al,  1993],  In  this  paper,  we  present  the  impact  that  high-energy  power  law  tails 
have  on  calculations  of  the  ionospheric  effects  of  precipitating  protons. 
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•  P.H.  Doherty,  D.N.  Anderson,  J.  Eicher  and  J.A.  Klobuchar,  ‘Total  Electron  Content  Over 
the  Pan  American  Longitudes:  March-April  1994”,  presented  at  the  Ionospheric  Effects 
Symposium  held  in  Alexandria,  VA  in  May,  1996.  It  also  appeared  in  1996  Ionospheric 
Effects  Symposium,  441-449, 1996. 
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Abstract 


An  experimental  campaign  to  measure  diurnal  changes  in  Total  Electron  Content  (TEC), 
over  the  wide  latitude  range  from  approximately  50  degrees  North  to  40  degrees  South  was 
carried  out  from  28  March  through  11  April  1994  by  monitoring  the  differential  carrier  phase 
from  the  US  Navy  TRANSIT  Navigation  Satellites,  using  a  chain  of  ground  stations  aligned 
along  the  approximate  70  degree  west  longitude  meridian.  This  campaign  was  conducted 
primarily  to  study  the  day-to-day  variability  of  the  equatorial  anomaly  region.  The  experimental 
plan  included  using  the  received  values  of  TEC  from  the  chain  of  stations  to  construct  profiles  of 
electron  density  versus  latitude  using  tomographic  reconstruction  techniques,  and  then,  to 
compare  these  reconstructions  against  a  theoretical  model  of  the  low-latitude  ionosphere. 

The  diurnal  changes  in  TEC  along  this  latitude  chain  of  stations  showed  a  high  degree  of 
variability  from  day-to-day,  especially  during  a  magnetic  storm  which  occurred  near  the 
beginning  of  the  campaign.  The  equatorial  anomaly  in  TEC,  showed  large  changes  in  character 
in  the  two  hemispheres,  as  well  as  differences  in  magnitude  from  day-to-day.  The  latitudinal 
gradients  of  TEC,  especially  in  the  lower  mid-latitudes,  also  showed  large  differences  between 
magnetic  storm  and  quiet  conditions.  Comparisons  of  the  TEC  data  with  the  theoretical  model 
illustrate  the  sensitivity  of  the  model  calculations  to  changes  in  magnetic  ExB  drift,  and  serves  to 
validate  the  strong  dependence  of  these  drifts  of  the  formation  and  the  strength  of  the  equatorial 
anomaly  regions. 

•  A.J.  Preble,  D.T.  Decker,  and  D.N.  Anderson,  “Improving  IRI90  Low  Latitude  Ionospheric 
Specification”,  presented  at  the  Ionospheric  Effects  Symposium  held  in  Alexandria,  VA  in 
May,  1996.  It  also  appeared  in  1996  Ionospheric  Effects  Symposium,  243-251, 1996. 

Abstract 

At  low  latitudes,  a  number  of  comparisons  between  the  IRI90  model  of  F-region  electron 
density  profiles  with  observed  profiles  measured  by  the  Jicamarca  incoherent  scatter  radar 
indicate  that  during  the  daytime,  the  observed  profile  shape  is  much  broader  in  altitude  than  that 
specified  by  IRI90;  while  at  night,  just  after  sunset,  observed  Hmax  values  are  significantly 
higher.  This  is  especially  true  during  periods  of  high  solar  activity.  The  theoretically-derived 
ionospheric  parameters  such  as  Hmax,  Nmax  and  profile  shape  which  are  contained  in  the 
Parameterized  Ionospheric  Model  (PIM)  have  been  shown  to  be  in  better  agreement  with 
Jicamarca  observations.  This  paper  describes  an  attempt  to  improve  IRI90  at  low  latitudes  by 
calculating  three  ionospheric  parameters,  the  profile  half  thickness,  Nmax,  and  Hmax  from  PIM 
as  a  function  of  latitude,  longitude,  and  solar  local  time  for  a  variety  of  solar  cycle,  seasonal  and 
geomagnetic  conditions.  The  generation  of  electron  density  profiles  using  these  three  parameters 
will  be  presented  as  well  as  a  description  of  how  these  expressions  might  be  implemented  into 
the  IRI90  model.  Finally,  we  discuss  the  degree  of  improvement  which  has  been  achieved  by 
comparing  these  profile  parameters  with  a  number  of  ionospheric  observations  not  only  at  the 
magnetic  equator  but  also  at  higher  magnetic  latitudes  near  the  crests  of  the  equatorial  anomaly. 


•  P.J.  Sultan,  “Satellite  Signatures  of  the  Global  Occurrence  Morphology  of  Equatorial  Spread- 
F\  presented  at  the  Ionospheric  Effects  Symposium  held  in  Alexandria,  VA  in  May,  1996.  It 
also  appeared  in  1996  Ionospheric  Effects  Symposium,  287-292, 1996. 
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•  P.J.  Gendron,  P.H.  Doherty,  and  J.A.  Klobuchar,  “Absolute  Real-Time  Ionospheric  Measure¬ 
ments  from  GPS  Satellites  in  the  Presence  of  Anti-Spoofing”,  presented  at  the  ION-GPS 
52nd  Annual  Meeting  held  in  Cambridge,  MA  in  June  1996.  It  also  appeared  in  the 
Proceedings  of  the  ION  52nd  Annual  Meeting,  547-556,  1996. 

Abstract 

Determining  user  position  and  velocity  using  signals  from  the  constellation  of  GPS 
satellites  depends  upon  receiving  signals  from  satellites  in  widely  disparate  directions  in  the  sky 
to  minimize  the  effects  of  geometric  dilution  of  precision  (GDOP).  The  range  error  due  to  the 
electron  content  of  the  earth’s  ionosphere  is  normally  corrected  by  using  signals  of  both  the  GPS 
L2  and  LI  frequencies  to  automatically  correct  for  the  first  order  ionospheric  range  error  effect. 
When  a  GPS  satellite  is  just  rising  above  the  horizon  it  is  desirable  to  use  it  in  the  navigation 
solution,  as  soon  as  it  reaches  as  elevation  above,  say  5°,  so  that  the  newly  acquired  satellite  aids 
in  producing  a  position  solution  with  a  minimum  GDOP.  However,  the  lower  signal-to-noise 
ratio  when  a  GPS  satellite  is  viewed  near  the  horizon,  especially  evident  when  anti-spoofing  (a/s) 
is  turned  on,  can  limit  the  accuracy  of  making  measurements  of  absolute  ionospheric  range  delay 
to  values  higher  than  approximately  one  meter  until  the  satellite  reaches  an  elevation  angle 
greater  than  approximately  30°. 

The  comparative  effects  of  a/s  and  high  multi-path  on  the  GPS  differential  group  delay 
signal  have  been  examined  and  the  results  clearly  demonstrate  that  even  in  a  high  multi-path 
environment,  a/s  is  the  dominant  effect  which  limits  early  determination  of  absolute  ionospheric 
range  error.  A  novel,  non-recursive,  robust  residual-adaptive  filter  was  found  to  give  the  best 
results  in  obtaining  absolute  ionospheric  range  error  corrections  at  low  elevation  angles  in  the 
presence  of  a/s.  The  times  required  from  satellite  rise  until  an  absolute  ionospheric  range  error  of 
less  than  0.16  meters  rms.  could  be  obtained  corresponded  to  elevation  angles  of  approximately 
30°  in  the  presence  of  a/s,  and  considerably  lower  with  a/s  turned  off.  Multi-path  was  a 
significant  limitation  in  determining  absolute  ionospheric  range  error  only  with  a/s  turned  off, 
and  then  only  in  a  relatively  bad  multi-path  environment. 


•  J.A.  Vladimer,  V.R.  Ewell,  M.C.  Lee,  P.H.  Doherty,  D.T.  Decker,  and  D.N.  Anderson, 
“Ionospheric  TEC  Observations  from  the  TOPEX  Satellite”,  presented  at  the  23rd  IEEE 
International  Conference  of  Plasma  Science  in  Boston,  MA  in  June  1996  and  the  PRIMO 
session  at  the  Cedar  Workshop  in  Boulder,  CO  in  June  1996. 

Abstract 

Variability  of  Total  Electron  Content  (TEC)  in  the  equatorial  anomaly  region  of  the 
ionosphere  can  be  studied  extensively  using  the  results  of  measurements  taken  by  the 
NASA/CNES  satellite,  TOPEX/Poseidon,  which  provides  global  ocean  coverage  (0  to  360° 
longitude,  -66  to  66°  latitude).  The  NASA  Radar  Altimeter  (NRA)  is  the  first  space-borne  dual- 
frequency  altimeter  capable  of  accurately  measuring  vertical  ionospheric  TEC  below  1340  km. 
TOPEX  TEC  observations  have  already  been  used  to  support  results  from  an  ionospheric 
measurements  campaign  that  was  conducted  in  the  equatorial  region  of  South  America  by 
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Phillips  Laboratory  in  the  spring  of  1994.  The  best  agreement  in  TEC  values  is  seen  during 
intervals  of  longitudinal  proximity  of  the  satellites’  paths. 

Vertical  measurements  from  TOPEX  are  obtained  at  a  rate  of  one  per  second  at  5.3  and 
13.6  GHz  along  the  same  surface  tracks  every  10  days  during  which  global  ocean  coverage  is 
accomplished.  TOPEX  data  records  from  9/92  to  9/95  have  been  processed  and  reduced  to  one 
CDROM  disc  using  one  degree  averaging  over  1 10  km  along  the  track  for  evaluation  of  the  TEC 
measurements.  Various  plotting  programs  have  been  developed  to  display  the  observations.  To 
facilitate  further  TEC  studies,  interface  software  has  been  prepared  to  provide  identification  of 
anomaly  formation  and  structure.  The  TOPEX  over-ocean  data  can  be  used  as  a  supplement  to 
land-based  measurements  in  applications  to  ionospheric  research  at  low  and  mid  latitudes.  This 
study  focuses  on  comparisons  between  TOPEX  vertical  TEC  data  and  GPS  equivalent  vertical 
TEC  measurements  taken  near  the  east  and  west  coastal  regions  of  South  America.  Also,  the 
Phillips  Laboratory  Global  Parameterized  Ionospheric  Model  (PIM)  is  utilized  in  an  effort  to 
estimate  slant-to-vertical  conversion  errors. 

•  P.J.  Sultan,  “DMSP-Based  Prediction  of  Equatorial  Scintillation”,  presented  at  the  Cedar 
Workshop  held  in  Boulder,  CO  in  June,  1996. 


Doherty,  D.T.  Decker,  D.N.  Anderson  and  B.D.  Wilson,  "Observed  Ionospheric 
Dependence  on  Solar  Activity:  Implications  for  a  New  Single  Frequency  GPS  User 
Algorithm  ,  presented  at  the  ION-GPS  96  Meeting  held  in  Kansas  City,  MO  in  September, 
1996.  It  also  appeared  in  the  Proceedings  of  The  9th  International  Technical  Meeting  of  The 
Satellite  Division  of  the  Institute  of  Navigation,  1996. 

Abstract 

The  ionosphere  is  the  most  substantial  source  of  error  in  GPS  positioning  and  navigation. 
Radio  waves  propagating  through  the  ionosphere  suffer  an  additional  time  delay  as  a  result  of 
their  encounter  with  the  total  number  of  free  electrons  (TEC)  along  the  path  from  the  satellite  to  a 
ground  station.  The  ionosphere  is  a  highly  variable  and  complex  region  that  is  a  function  of 
geographic  location,  local  time,  season,  geomagnetic  activity  and  solar  activity.  For  the  dual¬ 
frequency  GPS  user,  the  range  error  caused  by  the  ionosphere  is  corrected  for  by  using  the  signals 
of  both  the  GPS  LI  and  L2  frequencies.  For  the  single  frequency  user,  the  Ionospheric 
Correction  Algorithm  (ICA)  was  developed  to  correct  for  approximately  50  percent  rms  of  the 
error.  The  ICA  was  developed  in  1975  using  the  Bent  Model  [Llewellyn  and  Bent,  1973],  a  state- 
of-the-art  empirical  ionospheric  model  that  was  based  on  data  recorded  during  an  average  solar 
cycle.  Solar  cycle  maximum  values  since  that  period  have  been  much  higher.  As  a  result,  the 
ICA  does  not  adequately  represent  periods  of  extremely  high  solar  maximum  activity. 

The  aim  of  this  paper  is  to  aid  in  the  development  of  a  new  GPS  single  frequency  user 
algorithm  by  defining  an  improved  fit  between  long  term  solar  activity  and  to  determine  if  that  fit 
changes  for  different  seasons  and  geographic  locations.  This  study  is  accomplished  by 
correlating  ionospheric  measurements  recorded  at  several  widely  spaced  geographic  locations 
over  the  last  four  solar  cycles  with  various  proxies  of  the  10.7  cm  solar  flux  (Fi0.7).  Results  are 
provided  in  a  statistical  format  for  geomagnetically  quiet  and  disturbed  conditions.  The  results  of 
this  study  will  be  used  in  the  development  of  new  coefficients  that  describe  worldwide 
ionospheric  behavior.  The  improved  Ionospheric  Correction  Algorithm  will  benefit  all  users  of 
single-frequency  GPS  receivers. 
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•  C.G.  Stergis,  "Instrumentation  for  UV  and  Visible  Measurements  Above  Thunderstorms", 
presented  at  the  annual  meeting  of  the  International  Society  for  Optical  Engineering  held  in 
Denver,  CO  in  August,  1996.  Also  in  Proc.  Intern.  Soc.  Optical  Eng.,  Vol.  2831,  Nov.  1996. 


•  D.T.  Decker,  C.E.  Valladares,  and  D.N.  Anderson,  "Modeling  F-Region  Polar  Cap  Patches 
and  Boundary  Blobs",  presented  at  the  31st  COSPAR  Scientific  Assembly  held  in 
Birmingham,  United  Kingdom  in  July,  1996. 

Abstract 

Over  the  last  30  years,  observations  have  shown  that  the  high-latitude  F-region  contains  a 
variety  of  mesoscale  structures.  In  particular  by  the  end  of  the  1980s,  enough  was  known  that  a 
sequence  of  events  could  be  described  where  various  mesoscale  structures  evolved  from  one  to 
another.  The  sequence  begins  at  subauroral  latitudes  in  the  afternoon  sector  where  a  broad 
wedge  of  enhanced  ionization  is  observed.  At  higher  latitudes,  this  wedge  thins  in  longitude  to 
what  is  typically  referred  to  as  a  tongue  of  ionization.  Then  by  some  process,  that  wasn't  clearly 
defined  but  presumably  took  place  near  or  within  the  cusp  or  throat,  the  tongue  is  structured  into 
patches  which  then  drift  across  the  polar  cap.  The  patches,  in  turn,  are  then  reconfigured  into 
boundary  and  subauroral  blobs  as  the  enhanced  plasma  convects  out  of  the  polar  cap  into  the 
auroral  zone.  For  the  most  part,  this  picture  was  built  up  from  observations  with  little  theoretical 
support  from  F-region  modeling.  However,  over  the  last  several  years  there  have  been  significant 
successes  in  the  modeling  of  both  patches  and  blobs.  We  will  discuss  the  theoretical  approaches, 
simulations  involving  time-varying  global  convection,  simulations  involving  meso-scale 
convection  ( plasma  jet,  vortices,  flow  channels,  etc.),  and  future  directions. 


•  D.N.  Anderson,  D.T.  Decker,  and  P.H.  Doherty,  "Global  Ionospheric  Modeling  of 
Geomagnetic  Storms",  presented  at  the  31st  COSPAR  Scientific  Assembly  held  in 
Birmingham,  United  Kingdom  in  July,  1996. 

Abstract 

It  is  well  established  that  the  ionospheric  impact  of  geomagnetic  storms  is  global  in  nature; 
however,  many  past  theoretical  studies  have  focused  on  particular  sub  regions  of  the  ionosphere. 
For  example,  we  have  previously  used  a  low  latitude  ionospheric  model  to  study  the  impact  of 
storms  on  the  low-latitude  ionosphere.  In  particular,  we  used  as  input  to  the  model  the  observed 
vertical  drift  velocities  measured  by  the  Jicamarca,  Peru  incoherent  scatter  radar  facility  during 
storms  on  22,23  March  1990  and  13,14  June  1991.  The  resultant  simulations  produced  dramatic 
changes  in  the  electron  density  profiles  that  were  consistent  with  DMSP  in-situ  electron  density 
measurements.  This  same  model  has  also  been  used  in  other  studies  that  focused  on  specific 
storm  effects  in  the  low  and  mid-latitude  ionosphere.  However,  in  order  to  study  the  storm 
effects  on  a  global  scale,  it  would  be  useful  to  have  a  global  model.  In  recent  years,  we  have 
developed  such  a  global  theoretical  ionospheric  model  (GTIM)  and  successfully  used  it  to  study 
problems  at  high,  mid,  and  low  latitudes.  In  this  paper,  we  will  review  earlier  low  latitude  work 
as  well  as  discuss  our  initial  simulations  of  the  global  ionospheric  impact  of  geomagnetic  storms. 
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In  particular,  we  will  focus  on  the  uniqueness  of  individual  storms  and  the  difficulty  of  defining  a 
typical  or  average  storm. 


•  D.N.  Anderson,  D.T.  Decker,  and  C.E.  Valladares,  "Modeling  High  Latitude  F-reeion 
Mesoscale  Structures",  presented  at  the  25th  General  Assembly  of  URSI  held  in  Lille  France 
from  August  28  to  September  5, 1996. 


Abstract 

Over  the  last  30  years,  observations  have  shown  that  the  high-latitude  F-region  contains  a 
variety  of  mesoscale  structures.  In  particular,  it  has  been  found  that  during  periods  of  negative  or 
southward  directed  EL  enhanced  "patches"  of  ionization  drift  across  the  polar  cap  in  an 
antisunward  direction.  These  polar  cap  patches  come  in  a  wide  variety  of  shapes  and  scale  sizes 
The  typical  sizes  range  from  a  few  hundred  to  -1000  km,  and  the  enhancements  of  plasma 
densities  can  be  up  to  an  order  of  magnitude  above  the  surrounding  background.  Over  the  years, 
there  have  naturally  been  a  variety  of  suggestions  as  to  how  patches  might  be  produced;  however! 
identification  of  the  precise  mechanisms  of  patch  creation  has  remained  elusive. 

By  the  end  of  the  1980s,  enough  was  known  that  a  "descriptive  working  model"  could  be 
constructed.  This  model  described  a  sequence  of  events  where  various  large  scale  structures 
evolved  from  one  to  another.  The  sequence  begins  at  subauroral  latitudes  in  the  afternoon  sector 
where  a  broad  wedge  of  enhanced  ionization  is  observed.  At  higher  latitudes,  this  wedge  thins  in 
longitude  to  what  is  typically  referred  to  as  a  tongue  of  ionization.  Then,  by  some  process  that 
wasn't  clearly  defined,  but  presumably  took  place  near  or  within  the  cusp  or  throat,  the  tongue  is 
structured  into  patches  which  then  drift  across  the  polar  cap.  The  patches,  in  turn,  are  then 
reconfigured  into  boundary  and  subauroral  blobs  as  the  enhanced  plasma  convects  out  of  the 
polar  cap  into  the  auroral  zone. 

For  the  most  part,  this  picture  was  built  up  from  observations  with  some  theoretical  support 
from  F-region  modeling.  But  there  were  many  gaps  in  the  details,  many  uncertainties  and 
actually  very  few  theoretical  studies  had  explicitly  focused  on  how  F-region  structure  might  be 
produced.  However,  over  the  last  several  years,  modelers  have  successfully  simulated  the 
creation  of  high-latitude  large  scale  plasma  structure  using  two  separate  mechanisms  (time- 
varymg  global  convection  and  meso-scale  events).  Morphological  predictions  have  been  made 
that  are  currently  being  tested  against  data.  A  study  at  Sondrestrom  has  suggested  that  both 
mechanisms  can  be  important  and  may  very  well  operate  together  in  a  complex  and  intertwined 
manner.  Simulations  using  time-dependent  convection  have  also  illustrated  how  patches  can 
evolve  into  blobs  in  much  the  way  predicted  in  the  work  of  Robinson  et  al.m.  The  purpose  of 
this  paper  is  to  report  on  this  recent  progress  in  patch  and  blob  modeling.  We  will  discuss  the 
theoretical  approaches,  simulations  involving  time-varying  global  convection,  simulations 
involving  meso-scale  convection  (plasma  jet,  vortices,  flow  channels,  etc.),  and  future  directions. 

[1]  R.T.  Tsunoda,  Rev.  Geophys.,  26,  719, 1988. 

[2]  R.M.  Robinson  et  al.,  J.  Geophys.  Res.,  90, 7533, 1985. 


•  D.T.  Decker  and  D.N.  Anderson,  "Simulations  of  GPS/MET  Ionospheric  Observations", 
presented  at  the  25  General  Assembly  of  URSI  held  in  Lille,  France  from  August  28  to 
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September  5, 1996.  Also  presented  at  the  1996  Fall  AGU  meeting  held  in  San  Francisco,  CA 
in  December,  1996. 


Abstract 

GPS/MET  is  a  program  using  a  low-earth  orbiting  satellite  instrumented  with  a  GPS 
receiver  to  perform  limb  soundings  of  the  atmosphere  using  radio  occultation  of  GPS  satellites. 
The  primary  purpose  of  the  program  is  to  perform  remote  sensing  of  the  troposphere  and 
stratosphere.  However,  the  limb-viewing  geometry  does  provide  a  unique  opportunity  for 
ionospheric  research.  The  useful  measured  quantity  in  this  case  is  the  Total  Electron  Content 
(TEC)  between  the  low-earth-orbiting  satellite  and  a  particular  GPS  satellite.  In  contrast  to 
ground  based  GPS  measurements,  the  limb  viewing  geometry  can  provide  the  type  of  vertical 
information  that  is  poorly  measured  by  ground-based  GPS  measurements. 

In  this  paper,  the  Global  Theoretical  Ionospheric  Model  (GTIM)  is  used  to  simulate  the 
TEC  as  observed  by  both  ground-based  GPS  receivers  and  GPS/MET.  The  model  calculates  the 
atomic  oxygen  ion  density  as  a  function  of  altitude,  latitude,  and  local  time.  It  includes  the 
effects  of  production  of  ionization  by  solar  extreme  ultraviolet  radiation  and  electron 
precipitation;  loss  through  charge  exchange  with  molecular  nitrogen  and  oxygen;  and  transport 
by  diffusion,  neutral  winds,  and  ExB  convection  drifts.  The  simulations  are  designed  to  contrast 
the  signatures  of  various  features  of  the  ionosphere  such  as  the  Appleton  anomaly,  the  high 
latitude  trough,  boundary  and  auroral  blobs,  and  polar  cap  patches.  Of  particular  interest  is  the 
contrast  between  the  ground-based  signatures  and  GPS/MET  signatures.  We  also  contrast  the 
TEC  given  by  a  realistic  ionosphere  to  that  given  by  simple  approximations  of  the  ionosphere. 


•  P.H.  Doherty,  R.  Loh,  and  D.N.  Anderson,  “Statistics  of  the  Spatial  and  Temporal  Variations 
in  Ionospheric  Range  Delay”,  submitted  to  ION-GPS  97  to  be  held  in  Kansas  City,  MO  in 
September,  1997. 


Abstract 

In  this  study,  we  present  the  statistics  of  the  spatial  and  temporal  variations  in  ionospheric 
range  delay  recorded  at  several  locations  in  the  continental  US,  Canada,  Alaska  and  Hawaii. 
These  measurements,  recorded  from  1993  through  early  1997,  were  obtained  from  the 
International  GPS  Geodynamics  Service  (IGS)  network,  managed  by  the  Jet  Propulsion 
Laboratory.  The  statistics  of  spatial  variations  are  based  on  simultaneous  measurements  made 
along  different  lines  of  sight  at  each  location.  Statistics  of  the  temporal  ionospheric  variations 
are  based  on  the  highest  elevation  measurement  at  particular  times  of  day  at  each  individual 
receiver  location.  Since  these  measurements  cover  only  solar  moderate  to  minimum  conditions, 
we  augment  the  statistics  on  the  temporal  variations  with  a  historical  data  base  of  ionospheric 
measurements  recorded  at  similar  locations  over  a  full  solar  cycle.  The  historical  data  is  based 
on  Faraday  rotation  measurements  of  VHF  signals  from  geostationary  satellites  of  opportunity. 
Unfortunately,  the  sites  included  in  the  historical  data  are  too  widely  spaced  to  infer  spatial 
variations  in  the  next  solar  maximum. 

This  analysis  is  intended  to  provide  a  characterization  of  ionospheric  variations  that  will  be 
encountered  in  the  Federal  Aviation  Administration’s  (FAA)  operational  Wide  Area 
Augmentation  System  (WAAS).  The  WAAS  is  intended  to  support  en-route  and  precision 
approach  flight  operations.  An  important  factor  in  the  operational  system  will  be  the  accuracy  of 
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the  vertical  ionospheric  delay  estimates  that  will  be  provided  to  a  user  aircraft  in  a  satellite 
broadcast  message.  These  estimates  of  vertical  ionospheric  delay  will  be  based  on  slant  dual¬ 
frequency  GPS  measurements  of  ionospheric  range  delay  made  at  24  reference  stations  (21  in  the 
continental  US,  1  in  Hawaii,  1  in  Alaska  and  1  in  Puerto  Rico).  Historically,  the  ionosphere  has 
proved  to  be  extremely  variable  in  nature.  Its  behavior  is  a  function  of  geographic  location,  local 
time,  season,  geomagnetic  activity  and  solar  activity.  In  the  continental  US,  measurements  of 
ionospheric  range  delay  at  midday  can  range  from  less  than  2  meters  at  solar  minimum  to  nearly 
20  meters  at  solar  maximum.  In  Hawaii,  daytime  observations  of  range  delay  as  high  as  30 
meters  have  been  recorded  during  a  solar  maximum  period.  Fortunately,  these  large  variations  in 
the  ionosphere  develop  slowly  as  a  function  of  long  term  changes  in  solar  activity  and  are  not 
likely  to  cause  degraded  accuracy  in  the  vertical  ionospheric  delay  estimates  for  the  WAAS.  The 
most  critical  problems  of  ionospheric  estimation  for  the  WAAS  will  occur  during  periods  of 
ionospheric  scintillation  and  during  periods  of  major  geomagnetic  storm  activity  when  the  spatial 
and  temporal  gradients  in  the  ionosphere  can  be  significantly  higher  than  usual.  This  study  will 
quantify  the  magnitude  and  frequency  of  these  large  spatial  and  temporal  gradients  observed  with 
dual-frequency  GPS  from  single  locations  during  solar  moderate  to  minimum  periods.  It  will 
also  provide  a  historical  reference  to  the  temporal  gradients  observed  during  a  previous  solar 
maximum  period. 


•  D.T.  Decker,  C.E.  Valladares,  and  D.N.  Anderson,  “Specifying  the  High  Latitude  F  Region 
Weather  submitted  to  IAGA  meeting  to  be  held  in  Uppsala,  Sweden  in  August  1997. 

Abstract 

Observations  have  shown  that  much  of  the  high-latitude  F-region  weather  is  caused  by 
small  scale  (irregularities)  and  large  scale  electron  density  structures.  These  same  observations 
also  suggest  that  there  is  an  intimate  cause  and  effect  relation  between  the  large  scale  structures 
(>10  km)  and  the  smaller  scale  (<10  km)  irregularities.  Theoretical  work  on  modeling  large  scale 
structures  has  recently  demonstrated  that  we  can  simulate  the  formation  and  evolution  of  polar 
cap  patches  and  boundary  blobs,  two  types  of  F-region  structures  that  are  observed  coincident 
with  smaller  scale  irregularities.  However,  it  has  not  been  demonstrated  that  we  can  accurately 
model  such  structures  at  a  particular  place  and  time.  In  this  paper,  we  report  on  efforts  to  use  the 
Global  Theoretical  Ionospheric  Model  (GTIM)  to  determine  what  is  required  to  accurately 
specify  the  appearance  of  patches  and  blobs.  We  will  particularly  focus  on  the  sensitivity  of 
model  results  to  details  of  the  time-varying  convection  pattern  and  the  need  for  a  weather 
specification  of  the  high-latitude  electric  field. 


•  D.N.  Anderson,  D.T.  Decker,  P.H.  Doherty,  and  M.W.  Fox,  “Modeling  the  Low  Latitude  F 
Region  Weather”  submitted  to  1996  IAGA  meeting  to  be  held  in  Uppsala,  Sweden  in  August, 
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Abstract 


The  low-latitude  F-region  is  known  as  a  very  dynamic  region  where  dramatic  plasma 
irregularities  can  occur  over  a  large  range  of  scale  sizes  and  amplitudes.  However,  it  is  also  a 
region  where  the  ambient  F-region  plasma  has  great  day-to-day  variability.  In  large  part,  this 
variability  comes  from  the  day-to-day  behavior  of  the  F-region  electric  field,  in  particular,  the 
zonal  field  (vertical  drift).  This  is  so  because  the  vertical  drift  is  a  key  factor  in  determining  the 
basic  character  of  the  low-latitude  F-region.  Thus,  a  primary  requirement  for  modeling  the  low 
latitude  space  weather  will  be  the  ability  to  monitor  the  electric  field  and  use  it  as  a  weather  input 
to  an  ionospheric  model.  In  this  work,  we  are  examining  several  observational  sources  (radar, 
digisonde,  and  magnetometer)  for  the  day-to-day  behavior  of  the  vertical  drift  and  will  be  using 
the  resultant  drifts  as  inputs  to  the  Global  Theoretical  Ionospheric  Model  (GTIM).  The 
usefulness  of  these  sources  will  be  judged  by  comparisons  between  GTIM  results  and  various 
ionospheric  measurements  such  as  digisonde,  radar,  GPS,  and  TOPEX  observations. 


•  D.T.  Decker,  C.E.  Valladares,  E.  MacKenzie,  D.N.  Anderson,  S.  Basu,  and  Su.  Basu, 
“Modeling  the  Climatology  and  Weather  of  Blobs”,  submitted  to  1996  IAGA  meeting  to  be 
held  in  Uppsala,  Sweden  in  August,  1997. 

Abstract 

Over  the  last  30  years,  observations  have  shown  that  the  auroral  F  region  contains 
mesoscale  structures  known  as  blobs.  While  these  regions  of  enhanced  plasmas  have  been 
categorized  into  several  types  and  a  variety  of  sources  have  been  postulated,  most  work  has  been 
observational  and  only  a  few  theoretical  efforts  have  been  reported  [Robinson  et  al.,  1985; 
Anderson  et  al.,  1996].  In  the  recent  work  of  Anderson  et  al.  [1996],  the  Global  Theoretical 
Ionospheric  Model  (GTIM)  was  used  to  model  the  formation  of  boundary  blobs  using  time- 
varying  convection.  The  structures  produced  had  many  of  the  features  associated  with  boundary 
blobs.  However,  the  next  step  is  to  demonstrate  that  both  the  general  climatology  of  blobs  can  be 
modeled  and  that  blobs  in  a  particular  region  at  a  particular  time  can  be  modeled.  We  will  report 
on  efforts  to  compare  theoretical  blob  morphology  to  observed  scintillation  morphology  and  first 
efforts  to  reproduce  blobs  observed  on  a  specific  day  by  the  Chatanika  radar.  Discussion  will 
focus  on  the  global  observations  needed  to  provide  the  model  inputs  as  well  as  the  observations 
with  which  to  constrain  the  model  outputs. 

Anderson,  D.N.,  D.T.  Decker,  and  C.E.  Valladares,  "Modeling  Boundary  Blobs  Using  Time  Varying  Convection", 
Geophys.  Res.  Lett.,  23, 579-582, 1996. 

Robinson,  R.M.,  R.T.  Tsunoda,  and  J.F.  Vickrey,  "Sources  of  F  region  ionization  enhancements  in  the  nighttime 
auroral  zone,  J.  Geophys.  Res.,  90, 7533-7546, 1985. 
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Copies  of  reprints  of  the  following  two  papers  appear  in  this  section: 

C.E.  Valladares,  D.T.  Decker,  R.  Sheehan,  and  D.N.  Anderson,  "Modeling  the  formation 
of  polar  cap  patches  using  large  plasma  flows",  Radio  Science,  31,  573-593  May-June 
1996. 

•  D.T.  Decker,  B.V.  Kozelov,  B.  Basu,  J.R.  Jasperse,  and  V.E.  Ivanov,  "Collisional 
Degradation  of  the  Proton-H  Atom  Fluxes  in  the  Atmosphere:  A  Comparison  of 
Theoretical  Techniques",  J.  Geophys.  Res.,  101,  26947-26960, 1996. 


A  copy  of  the  following  camera-ready  paper  also  appears  in  this  section: 

•  P.H.  Doherty,  D.N.  Anderson,  and  J.A.  Klobuchar,  “Total  Electron  Content  Over  the  Pan 
American  Longitudes:  March-April  1994”,  Radio  Science,  in  press. 
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Radio  Science ,  Volume  31,  Number  3,  Pages  573-593,  May-June  1996 


Modeling  the  formation  of  polar  cap  patches  using 
large  plasma  flows 

C.  E.  Valladares,  D.  T.  Decker,  and  R.  Sheehan 

Institute  for  Space  Research,  Boston  College,  Newton  Center,  Massachusetts 


D.  N.  Anderson 

Phillips  Laboratory,  Geophysics  Directorate,  Hanscom  Air  Porce  Base,  Massachusetss 


Abstract.  Recent  measurements  made  with  the  Sondrestrom  incoherent  scatter  radar  have 
indicated  that  the  formation  of  polar  cap  patches  can  be  closely  associated  with  the  flow  of  a 
large  plasma  jet.  In  this  paper,  we  report  the  results  of  a  numerical  study  to  investigate  the  role 
of  plasma  jets  on  patch  formation,  to  determine  the  temporal  evolution  of  the  density  structure, 
and  to  assess  the  importance  of  0+  loss  rate  and  transport  mechanisms.  We  have  used  a  time- 
dependent  model  of  the  high-latitude  F  region  ionosphere  and  model  inputs  guided  by  data 
collected  by  radar  and  ground-based  magnetometers.  We  have  studied  several  different 
scenarios  of  patch  formation.  Rather  than  mix  the  effects  of  a  complex  of  variations  that  could 
occur  during  a  transient  event,  we  limit  ourselves  here  to  simulations  of  three  types  to  focus  on  a 
few  key  elements.  The  first  attempt  employed  a  Heelis-type  pattern  to  represent  the  global 
convection  and  two  stationary  vortices  to  characterize  the  localized  velocity  structure.  No 
discrete  isolated  patches  were  evident  in  this  simulation.  The  second  modeling  study  allowed 
the  vortices  to  travel  according  to  the  background  convection.  Discrete  density  patches  were 
seen  in  the  polar  cap  for  this  case.  The  third  case  involved  the  use  of  a  Heppner  and  Maynard 
pattern  of  polar  cap  potential.  Like  the  second  case,  patches  were  seen  only  when  traveling 
vortices  were  used  in  the  simulation.  The  shapes  of  the  patches  in  the  two  cases  of  moving 
vortices  were  defined  by  the  geometrical  aspect  of  the  vortices,  i.e.  elliptical  vortices  generated 
elongated  patches.  When  we  “artificially”  removed  the  Joule  frictional  heating,  and  hence  any 
enhanced  0+  loss  rate,  it  was  found  that  transport  of  low  density  plasma  from  earlier  local  times 
can  contribute  to  -60%  of  the  depletion.  We  also  found  that  patches  can  be  created  only  when 
the  vortices  are  located  in  a  narrow  local  time  sector,  between  1 000  and  1 200  LT  and  at  latitudes 
close  to  the  tongue  of  ionization. 


1.  Introduction 

Large-scale  density  structures  are  commonly  observed 
in  the  polar  cap.  When  the  interplanetary  magnetic  field 
(IMF)  is  directed  southward  mesoscale  (100  to  1000 
km),  density  enhancements,  named  polar  cap  “patches”, 
drift  across  the  polar  cap  in  the  antisunward  direction 
[Buchau  et  al ,  1983;  Weber  et  aL ,  1984,  1986;  Fukui  et 
al .,  1994].  When  the  IMF  is  directed  northward, 
elongated  streaks  of  precipitation-  enhanced  F  region 
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plasma  extend  aligned  to  the  noon-midnight  meridian  and 
move  toward  dawn  or  dusk  [Buchau  et  al. ,  1983;  Carlson 
et  al. ,  1984;  Valladares  et  ai,  1994a].  Both  types  of 
structures  can  be  associated  with  intense  levels  of 
scintillation  [Weber  et  al 1984;  Buchau  et  al. ,  1985; 
Basu  et  al. ,  1985,  1989]  and  disrupt  radio  satellite 
communication  systems.  In  this  paper,  we  investigate 
theoretically  the  formation  of  plasma  density 
enhancements  that  occur  during  B,  south  conditions.  To 
conduct  this  study,  we  have  used  a  three-dimensional 
time-dependent  model  of  the  high-latitude  F  region 
ionosphere  developed  by  Anderson  et  ai  [1988,1996]  and 
Decker  et  ai  [1994],  Differing  from  previous  attempts, 
we  have  incorporated  into  the  model  analytical 
expressions  of  localized  electric  field  structures  and  their 
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time-dependence.  We  checked  the  validity  of  the  model 
results  by  conducting  a  detailed  comparison  with  data 
gathered  from  the  Sondrestrom  incoherent  scatter  radar 
(ISR)  during  patch  formation  events. 

The  source  of  the  plasma  and  the  physical  processes 
that  form  the  polar  cap  patches  have  been  under 
investigation  for  over  10  years.  Buchan  et  al.  [1985]  and 
dc  la  Beaujardiere  el  al.  [1985]  were  the  first  researchers 
to  address  the  question  of  the  origin  of  the  density  inside 
the  patches.  Buchau  et  al.  [1985]  measured  values  of/,F2 
at  Thule,  Greenland  (86°  magnetic  latitude),  showing 
large  fluctuations.  The  minimum  f,F2  values  were  equal 
to  densities  produced  locally  by  the  sun  EUV  radiation, 
but  the  maximum  values  were  similar  to  the  densities  that 
were  produced  at  locations  equatorward  of  the  auroral 
oval.  These  observations  suggested  that  the  plasma 
density  inside  the  patches  was  produced  by  solar  radiation 
in  the  sunlit  ionosphere,  probably  at  subcusp  latitudes, 
and  then  carried  into  and  across  the  polar  cap  by  the 
global  convection  pattern.  In  fact,  Weber  et  al.  [1984] 
and  Foster  and  Doupnik  [1984]  confirmed  the  presence 
of  a  large  eastward  electric  field  near  midday,  likely 
directing  the  subauroral  plasma  poleward.  Buchau  et  al. 
[1985]  found  also  that  the  occurrence  of  patches 
displayed  a  strong  diurnal  UT  control.  The  patches  were 
seen  at  Thule  almost  exclusively  between  1200  and  0000 
UT.  This  UT  control  of  the  polar  cap  F  layer  ionization 
was  interpreted  in  terms  of  the  displacement  of  the 
convection  pattern  with  respect  to  the  geographic  pole.  It 
is  because  of  this  displacement  that  the  polar  convection 
is  able  to  embrace  higher-density  plasma  only  at  certain 
UT  hours  [Sojka  et  al.,  1993,  1994],  However,  the 
mechanism  (or  mechanisms)  by  which  the  auroral  plasma 
can  break  up  into  discrete  entities  is  still  a  matter  of 
debate. 

Tsunoda  [1988]  summarized  the  role  of  different 
mechanisms  conducive  for  patch  generation.  He 
suggested  that  changes  in  the  By  and/or  B.  components  of 
the  IMF  could  originate  temporal  variations  in  the  global 
flow  pattern,  drastically  disturbing  the  density  distribution 
within  the  polar  cap.  In  the  mid  1980s,  several  theoretical 
studies  were  designed  to  study  the  production,  lifetime, 
and  decay  of  large-scale  structures  inside  the  polar 
ionosphere  [Sojka  and  Schunk,  1986;  Schunk  and  Sojka, 
1987],  These  authors  indicated  that  hard  precipitation 
could  produce  plasma  enhancements  similar  to  the 
patches  and  blobs  found  in  the  auroral  oval.  They  argued 
that  in  the  absence  of  solar  radiation,  although  the 
particle-produced  E  region  will  rapidly  recombine,  the 
longer  lifetime  enhanced  F  region  ionization  could  build 
up  and  persist  for  several  hours.  Sojka  and  Schunk 


[1988]  theoretically  demonstrated  that  large  electric 
fields  could  create  regions  of  density  depleted  by  a  factor 
of  4.  Anderson  et  al.  [1988]  presented  a  model  of  the 
high-latitude  F  region  that  was  used  to  investigate  the 
effect  of  sudden  changes  in  the  size  of  the  polar  cap  upon 
density  enhancements  transiting  the  polar  cap.  While 
these  model  studies  were  able  to  form  structures,  most  of 
the  success  in  patch  modeling  has  been  attained  only  after 
the  first  High  Latitude  Plasma  Structures  meeting  was 
convened  at  Peaceful  Valley,  Colorado,  in  June  1992. 
Sojka  et  al.  [1993]  conducted  numerical  simulations  of 
the  effect  of  temporal  changes  of  the  global  pattern.  The 
latter  two  studies  were  successful  in  producing  density 
structures  at  polar  cap  latitudes.  Decker  et  al.  [1994]  used 
six  different  global  convection  patterns  and  localized 
velocity  structures  to  reproduce  the  digisonde 
measurements  of f,F,  values  at  Sondrestrom.  Lockwood 
and  Carlson  [1992]  used  the  formulation  of  transient 
reconnection  [Cowley  et  al.,  1991]  to  suggest  that  a 
transient  burst  of  reconnection  together  with  the 
equatorward  motion  of  the  ionospheric  projection  of  the 
X  line  could  extract  a  region  of  the  high  density 
subauroral  plasma,  divert  the  subauroral  density 
poleward,  and  finally  "pinch  off'  the  newly  formed  patch. 
More  recently,  Valladares  et  al.  [  1 994b]  has  shown  a  case 
study  where  a  fast  plasma  jet  containing  eastward  directed 
velocities  in  excess  of  2.5  km  s  1  was  able  to  increase  the 
O  recombination  rate  and  yield  an  east-west  aligned 
region  of  reduced  densities  across  a  poleward  moving 
tongue  of  ionization  (TOI).  Rodger  et  al.  [1994] 
presented  data  from  the  Halley  Polar  Anglo-American 
Conjugate  Experiment  radar  in  Antarctica,  suggesting  that 
a  region  of  depleted  densities  can  be  carved  out  by 
increased  0+  recombination  due  to  large  plasma  jets. 

This  paper  presents  a  detailed  modeling  of  the 
mechanism  of  patch  formation  presented  by  Valladares  et 
al.  [1994b]  and  Roger  et  al.  [1994].  The  mechanism 
described  here  operates  in  conjunction  with  a  background 
poleward  convection  to  produce  mesoscale  polar 
structures  from  subauroral  plasma.  This  process  may  also 
have  an  essential  role  in  determining  the  size  and  the 
shape  of  the  polar  cap  patches.  The  basic  elements  of  this 
mechanism  are  (1)  a  fast  plasma  jet,  observed  near  the 
midday  auroral  oval,  (2)  enhanced  ion  temperature  due  to 
Joule  frictional  heating,  (3)  enhanced  recombination  rates 
of  0+,  and,  (4)  transit  of  low  density  plasma  from  earlier 
local  times. 

The  paper  has  been  organized  in  the  following  manner. 
Section  2  succinctly  describes  our  model  of  the  high- 
latitude  ionosphere  [Anderson  el  al.,  1988;  Decker  et  al., 
1994]  and  the  calculations  of  the  localized  electric  fields 
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incorporated  into  the  model.  Section  3  presents  results  of 
a  computer  model  of  a  patch  formation  mechanism  using 
convection  patterns  for  typical  IMF  By  positive  and 
negative  values.  An  assessment  of  the  ability  of 
enhanced  recombination  loss  to  erode  parts  of  the  TOI 
and  to  disconnect  regions  from  the  oval  is  also  discussed 
in  this  section.  Section  4  presents  model  calculations  for 
vortices  located  at  several  different  local  times.  The 
paper  concludes  with  the  discussion  and  conclusions 
section. 

2.  Model  Description 

The  Global  Theoretical  Ionospheric  Model  (GTIM) 
calculates  the  density  altitude  profile  of  a  single  ion 
(0+)  along  a  flux  tube.  It  solves  the  coupled  continuity 
and  momentum  equations  for  ions  and  electrons.  The 
solution  of  the  differential  equations  is  simplified  by 
selecting  a  coordinate  system  in  which  one  dimension  is 
defined  parallel  to  the  local  magnetic  field.  The  model 
attains  three-dimensionality  by  repeating  the  profile 
calculations  along  many  (a  few  thousand)  flux  tubes. 
The  locations  of  the  flux  tubes  are  selected  to  cover  the 
range  of  latitudes  and  local  times  desired  for  the 
simulation.  The  mathematical  foundation  of  the  GTIM 
model  was  introduced  by  Anderson  [1971,1973]. 
Although  initially  used  to  study  the  low-latitude 
ionosphere,  it  has  been  extended  to  include  physical 
processes  of  the  high-latitude  ionosphere,  such  as  the 
effects  of  large  electric  fields  and  particle  precipitation 
[Anderson  et  al. ,  1988;  Decker  et  al ,  1994].  Recently, 
Decker  et  al  [1994]  have  described  how  the  inputs  to 
the  GTIM  model  can  easily  accommodate  a  time- 
dependent  convection  pattern  and  spatially  localized 
regions  containing  high  flows  to  represent  highly 
variable  situations,  such  as  those  existing  in  the  cusp 
region.  Other  geophysical  input  parameters,  such  as  the 
neutral  density  and  wind,  along  with  several  initial 
conditions,  such  as  the  plasma  density,  can  be  freely 
defined  to  simulate  different  scenarios.  In  this  paper, 
we  have  adjusted  the  global  convection  pattern  and 
included  a  local  electric  field  to  reproduce  the  velocities 
measured  on  February  19,  1990,  by  the  Sondrestrom 
ISR. 

The  radar  and  magnetometer  measurements  used  to 
define  the  model  inputs  were  discussed  by  Valladares  et 
al  [1994b].  These  authors  observed  a  large  channel 
oriented  in  the  east-west  direction  containing  jet-type 
eastward  velocities  of  order  2  km  s  \  In  this  channel, 
the  ion  temperature  (7))  was  enhanced  and  the  density 
(Nc)  was  depleted.  These  signatures  in  the  Tf  and  Ne 


values  implied  a  likely  increase  in  the  O  recombination 
rate.  Successive  radar  scans  indicated  that  the  plasma 
jet  was  continuously  moving  poleward.  This  motion 
was  also  confirmed  by  the  large  negative  deflection 
seen  at  later  times  by  magnetometers  located  at  latitudes 
poleward  of  Sondrestrom.  The  large  negative  deviation 
of  the  H  component  due  to  Hall  currents  appeared  at 
Qaanaaq  29  min  after  being  observed  at  Sondrestrom. 
This  implied  an  average  poleward  displacement  of  620 
m  s*1  for  the  plasma  jet.  Valladares  et  al  [1994b]  also 
presented  equivalent  velocity  vectors  deduced  from  the 
magnetometers  located  along  the  east  coast  of 
Greenland,  implying  a  flow  vorticity.  Similar  vorticity 
was  seen  in  the  resolved  radar  velocity  vectors.  The 
radar  line-of-sight  velocity  also  indicated  that  adjacent 
to  and  both  northward  and  southward  of  the  large 
plasma  jet  there  existed  regions  of  westward  flows.  In 
summary,  the  radar  and  magnetometer  data  suggested 
the  presence  of  two  adjacent  vortices  of  opposite 
vorticity  with  the  common  region  in  the  middle 
comprising  the  plasma  jet.  The  GTIM  model  also 
requires  other  geophysical  parameters,  such  as  the 
neutral  density,  wind  and  temperature,  the  ionization, 
and  chemical  loss  rates,  and  the  diffusion  coefficients. 
These  latter  atmospheric  parameters  were  selected  as 
described  by  Decker  et  al  [1994].  The  neutral  densities 
and  winds  were  calculated  using  the  mass 
spectrometer/incoherent  scatter  1986  (MSIS-86)  model 
[. Hedin ,  1987]  and  Hedin’s  wind  model  [Hedin  et  al , 
1991].  The  ion  loss  rate  was  computed  as  a  function  of 
an  effective  temperature.  This  parameter  is  derived 
using  the  following  equation  of  Schunk  et  al  [1975]: 

Teff  =  Tn  +  0.329E2 

where  E  is  the  magnetospheric  electric  field  in 
millivolts  per  meter. 

Plate  1  displays  the  NmF2  values  of  the  high  latitude 
ionosphere  as  a  function  of  magnetic  latitude  and  local 
time.  These  values  were  obtained  with  the  GTIM 
model  after  following  7200  flux  tubes  during  8  hours  of 
simulation  time.  During  this  time,  the  global 
convection  pattern  and  other  relevant  ionospheric 
parameters  were  kept  constant.  The  purpose  here  was 
to  obtain  the  initial  ionospheric  densities  to  be  used  as  a 
convenient  starting  point  for  simulating  the  ionospheric 
effects  of  a  plasma  jet.  The  actual  ionosphere  rarely  has 
8  hours  of  such  steady  conditions.  However,  this  allows 
us  to  single  out  individual  effects  that  can  be  produced 
by  various  ionospheric  processes.  The  real  ionosphere 
may  be  a  superposition  of  several  of  these  plasma  jet 
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Plate  I.  Polar  plot  of  the  NmF2  (peak  F  region  density)  of  the  high-latitude  ionosphere  at  1300  UT.  This 
density  plot  corresponds  to  the  initial  values  used  in  the  simulations  of  sections  3.1,  3.2,  3.4,  and  4.  The  values 
in  this  plot  were  obtained  by  running  the  model  for  8  hours  and  using  steady  state  inputs.  We  include  in  the 
lower  right  comer  the  black  and  white  version  of  this  plot  in  a  format  that  will  be  shown  in  subsequent  plots. 
The  black  dots  correspond  to  the  locations  of  the  Sondrestrom  and  Qaanaaq  stations,  and  the  white  dot  indicates 
the  location  of  the  European  Incoherent  Scatter  station.  The  latitudinal  circles  are  in  steps  of  10°. 


and  density  break-off  events.  A  prominent  feature  of 
Plate  1  is  the  presence  of  the  TOI.  At  1300  UT  it 
extends  from  longitudes  close  to  European  Incoherent 
Scatter  EISCAT  (19°E)  up  to  51°W,  where  it  turns 
poleward.  The  TOI  is  bounded  at  the  equatorward  edge 
by  a  region  of  densities  reduced  by  30%  and  poleward 
by  densities  almost  an  order  of  magnitude  smaller.  The 
larger  densities  in  the  TOI  are  due  mainly  to  two 
factors,  a  small  westward  flow  in  the  dusk  cell  and  an 
upward  lift  of  the  F  region.  The  longer  transit  time  and 
the  relatively  smaller  solar  zenith  angle  permit  the  solar 
radiation  to  build  up  much  higher  densities  in  this 


confined  region.  Similar  steady  structures  have  been 
presented  by  Crain  et  al.  [1993]  in  their  total  electron 
content  maps  of  the  high-latitude  region. 

The  vorticity  suggested  by  the  radar  and 
magnetometer  data  and  other  theoretical  implications  of 
solar  wind  -  magnetosphere  interactions  [Newell  and 
Sibeck,  1993]  led  us  to  infer  that  the  plasma  jet  may 
actually  consist  of  a  system  of  two  vortices 
superimposed  on  the  background  convection.  With  this 
as  a  guide,  we  searched  for  a  system  of  two  ellipsoidal 
potential  vortices  added  to  a  global  convection  pattern. 
The  search  was  carried  out  by  varying  the  cross  polar 
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Table  1.  Global  and  Local  Velocity  Patterns 


Polar  Cap 

Number  of 

Vortex 

Latitude  of 

Local  Time  of 

Global  pattern 

By  Potential, 
kV 

Vortices 

Major  Axis, 
deg 

Minor  axis, 
deg 

Potential, 

kV 

Vortices, 

deg 

Vortices, 

deg 

Heelis  12°  radius 

+  80 

2 

10,  10 

2.4,  1 

20,  -5 

74,  71.5 

11,  11 

Heppner-Maynard 

76 

2 

10,  10 

2.4,  1 

20,  -5 

74,  71.5 

11,  11 

cap  potential,  the  global  pattern,  the  radius  of  the  polar  convection  pattern.  Figures  la  and  lb  show  the  Heelis- 
cap  and  other  geometrical  parameters  of  the  two-vortex  type  global  convection  pattern  and  the  system  of  the 

system.  The  size,  peak-to-peak  voltage,  location,  two  vortices  of  Table  1.  Figure  lc  depicts  the  result  of 

orientation  and  aspect  of  the  vortices  were  adding  the  velocity  patterns  of  Figures  la  and  lb. 

systematically  changed  iteratively  to  obtain  the  best  fit  Figure  Id  represents  the  simulated  line-of-sight 

to  the  radar  velocities.  We  tried  over  10  million  velocities  which  would  have  been  measured  at 
combinations.  Table  1  presents  the  general  parameters  Sondrestrom  if  the  pattern  in  Figure  lc  had  been  in 

of  the  vortices  and  the  global  pattern  that  provided  the  effect.  The  agreement  between  Figure  Id  and  Figure  6a 

best  fit  to  the  Sondrestrom  velocity  data.  A  Heelis-type  of  Valladares  et  al.  [1994b],  reproduced  here  as  Figure 

pattern  gave  the  smallest  error.  However,  a  Heppner  le,  is  very  good.  Figure  If  shows  the  total  vector 

and  Maynard  By< 0  pattern  provided  errors  almost  as  velocity  that  was  obtained  using  the  convection  pattern 

small  as  the  Heelis-type  pattern.  Because  the  radar  of  Figure  lc. 

azimuthal  scan  covers  a  small  region  of  the  total  pattern,  Figure  2  shows  the  plasma  density  measured  by  the 

this  method  is  not  very  sensitive  in  discriminating  the  Sondrestrom  ISR  on  February  19,  1990,  during  four 


Figure  1.  Series  of  polar  plots,  (a)  High-latitude  potential  corresponding  to  a  Heelis-type  pattern  for  By  «  +8 
nT.  The  cross  polar  cap  potential  is  80  kV,  and  the  radius  of  the  polar  cap  is  12°.  (b)  the  potential  of  the 
convective  vortices,  (c)  Addition  of  the  potentials  in  Fgures  la  and  lb.  (d)  Simulated  Sondrestrom  line-of-sight 
velocities  considering  that  the  potentials  of  Figure  lc  were  in  effect,  (e)  Line-of-sight  velocities  measured  by  the 
Sondrestrom  radar  on  February  19,  1990  [Valladares  et  al.,  1994b,  Figure  6a].  (f)  Simulated  vector  velocities 
using  the  potential  of  Figure  lc. 
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West  Distance  (km)  East 


consecutive  elevation  scans.  Figure  2a,  obtained 
between  1251  and  1254  UT  and  before  the  appearance 
of  the  plasma  jet,  shows  a  density  value  above  3x10s 
cm  to  the  south  of  the  radar.  Overhead  and  to  the 
north,  the  density  is  structured  and  peaks  at  5.5x10s  cm' 
.  These  values  are  in  good  qualitative  agreement  with 
the  modeled  densities  of  Plate  1.  The  maximum 
modeled  peak  densities  are  5x10s  cm'3  and  9x10s  cm'3 
for  the  south  and  north  locations  respectively. 


3.  Patch  Modeling  Using  Mesoscale 
Velocity  Structures 

During  the  last  decade,  a  large  array  of  transient  events 
have  been  measured  in  the  high-  latitude  ionosphere. 
Convective  vortices,  double  vortices,  east-west  elongated 
ellipses,  and  vortices  moving  poleward,  equatorward, 
dawnward  and  duskward  have  all  been  observed.  While 
the  magnetic  response  of  these  events  have  been  studied 
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in  some  detail  [Gocrtz  el  a!.,  1985;  Heikkila  el  a /.,  1989; 
Ma  et  at.,  1991],  very  little  is  known  about  their  effect 
upon  the  ionospheric  densities.  Furthermore,  transient 
events  have  been  observed  to  have  a  fairly  complex 
temporal  behavior  that  is  not  well  understood.  As  a 
result,  rather  than  attempting  a  simulation  that  tries  to 
include  all  of  the  observed  behavior  of  a  transient  event, 
we  have  performed  four  numerical  simulations  that  focus 
on  just  a  few  key  elements  of  the  event.  The  four 
scenarios  modeled  were  chosen  to  help  us  learn  how  the 
density  evolves  under  different  conditions  of  the  local  and 
global  patterns  and  to  give  us  some  insight  as  to  how 
sensitive  our  results  are  to  the  various  properties  of  these 
events.  In  the  next  four  subsections,  we  report  on  the 
results  from  these  simulations.  First,  we  describe  a  case 
study  in  which  a  pair  of  transient  vortices  are  kept  fixed 
in  a  corrected  geomagnetic  coordinate  system.  In  the 
next  case,  we  include  the  observation  that  these  transient 
vortices  move  by  allowing  the  vortices  to  drift  according 
to  the  background  polar  cap  flow.  We  continue  with  a 
model  case  study  that  uses  a  Heppner-Maynard-type 
convection  pattern  as  the  global  convection.  This  allows 
us  to  assess  the  sensitivity  of  our  previous  results  to  the 
global  pattern  that  is  chosen.  Finally,  in  the  fourth 
simulation,  we  establish  the  importance  of  the 
enhancement  of  the  recombination  rate  in  producing 
regions  of  depleted  density  across  the  TOl. 

3.1.  Stationary  Vortices 

The  first  modeling  effort  was  implemented  using  two 
vortices  which  were  maintained  stationary  in  a  local  time 
versus  latitude  frame  of  reference.  The  vortices'  potential 
is  set  equal  to  the  predetermined  value  at  1300  UT;  15 
min  later  the  potential  is  turned  off.  Figure  3  and 
successive  figures  present  a  limited  region  of  the  northern 
hemisphere  high  latitude  ionosphere.  Here  we  show  a 
series  of  snapshots  of  the  ArmF2  corresponding  to  nine 
different  times  during  the  formation  of  the  density 
structure.  The  area  covered  in  each  panel  is  about  3500  x 
3500  km2.  Figures  3a  through  3f  are  3  min  apart.  Figures 
3g  -  3i  display  the  density  in  steps  of  10  min.  Figure  3a 
displays  the  initial  density  at  the  time  of  the  potential 
vortices  onset.  This  plot  corresponds  to  a  view  looking 
down  from  a  satellite  that  is  traveling  from  the  pole 
toward  the  equator.  Lower  latitudes  are  at  the  top  of  the 
figure,  and  poleward  latitudes  are  at  the  bottom.  The 
magnetic  pole  is  at  the  center  of  the  bottom  edge  of  each 
panel.  The  locations  of  the  Sondrestrom  and  Qaanaaq 
facilities  are  indicated  by  white  dots. 


The  TOI  is  clearly  depicted  in  Figure  3a;  it  extends 
nearly  aligned  with  the  equipotential  lines  and  turns 
poleward  near  Sondrestrom.  We  remind  the  reader  that 
the  N„,F2  pattern  of  Figure  3a  was  calculated  for  February 
19  at  1300  UT  using  the  actual  magnetic  and  solar  flux 
conditions  measured  on  that  day  and  the  inferred  size  of 
the  polar  cap,  as  defined  in  section  2.  A  weaker  TO!  may 
be  obtained  for  other  days  of  the  year,  different  UT,  or 
even  different  parameters  of  the  convection  pattern. 
Figure  3b  displays  the  peak  density  3  min  after  the  vortex 
potential  has  been  activated.  Even  at  this  time,  the  TOI  is 
already  drastically  distorted.  A  section  of  the  TOI  has 
rotated  anticlockwise;  simultaneously,  three  parallel 
channels  containing  depleted  densities  have  been  created, 
two  equatorward  and  one  poleward  of  the  fragmented 
TOI  segment.  All  three  elongated  traces  of  low  density 
are  coincidental  with  the  equatorward,  center,  and 
poleward  walls  of  the  vortices,  which  is  where  the  high 
speed  flows  and  the  enhanced  temperatures  are.  Three 
minutes  later  at  1306  (Figure  3c),  all  three  low-density 
channels  are  further  elongated  in  the  east-west  direction. 
A  region  of  low  density  starts  propagating  poleward  at  the 
poleward  edge  of  the  vortex  system.  The  following  three 
panels  (Figures  3d  to  3f)  show  further  elongation  and 
deeper  depletions  in  all  three  density  channels.  Clearly, 
the  presence  of  the  vortices  has  produced  density 
structures.  However,  the  density  structure  aligned  in  the 
east-west  direction  (Figure  3f)  is  still  attached  to  what 
remains  of  the  TOI.  After  the  localized  potential  is  turned 
off  at  1315  UT,  the  density  structures  lose  their  east- west 
alignment.  Figures  3g  to  3i  show  that  the  density 
enhancement  acquires  instead  a  more  north-south 
alignment.  This  additional  change  in  the  shape  of  the 
density  structure  is  due  to  a  nonuniform  velocity  near  the 
center  of  the  polar  cap  of  the  Heelis-type  pattern  in  Figure 
la.  If  the  flow  were  uniform,  then  the  equipotential  lines 
should  have  been  parallel  lines  in  the  polar  plot. 

Figure  4  presents  simulated  radar  elevation  scans 
through  the  modeled  ionospheric  volume.  Each  of  the 
panels  shows  an  equivalent  elevation  scan  in  the  meridian 
plane  that  would  have  been  measured  at  Sondrestrom  if 
the  modeled  densities  were  the  real  densities  around  the 
station.  This  scheme  allows  us  to  perform  a  one-to-one 
comparison  with  the  data  presented  by  Valladares  et  al. 
[1994b].  It  help  us  to  better  understand  the  radar  data 
when  only  elevation  scans  are  conducted  during  the 
experiment.  Each  of  the  scans  of  Figure  4  are  constructed 
from  the  density  data  displayed  in  the  corresponding 
panels  of  Figure  3.  The  orientation  of  the  scans  is  shown 
by  the  white  lines  in  each  panel  of  Figure  3.  The  larger 
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Figure  3.  Each  panel  presents  NmF2  values  of  a  section  of  the  high-latitude  ionosphere  at  different  times  during 
the  modeling  that  uses  stationary  vortices  as  input.  See  section  3.1  for  more  details.  The  white  dots  correspond  to 
the  locations  of  the  Sondrestrom  and  Qaanaaq  stations.  The  white  line  that  crosses  the  Sondrestrom  site 
indicates  the  extension  of  the  elevation  scans  of  Figure  2. 


density  at  the  northern  part  of  the  scan  of  Figure  4a  is  due 
to  the  passage  of  the  TOI.  Figure  4b,  corresponding  to  3 
min  along  the  simulation,  shows  mesoscale  structuring 
coincident  with  the  vortices.  Two  deep  depletions  are 
evident,  one  overhead  and  the  other  500  km  north  of  the 
station.  The  density  of  the  northern  depletion  is  3x10  cm* 
\  The  density  remains  unperturbed  at  the  northern  and 
southern  extremes  of  the  scan;  no  vortices  are  located 


there.  At  1309  UT,  9  min  after  the  beginning  of  the 
simulation,  both  depletions  reach  a  minimum  value  of 
2xl05  cm'3.  The  Ne  depletion  located  almost  overhead 
coincides  with  the  center  of  the  vortex,  where  the  plasma 
flow  is  maximum.  Figure  4e  presents  the  modeled 
densities  12  min  after  the  onset  of  the  potential  vortices. 
These  contours  can  be  closely  compared  with  Figure  2c 
which  shows  the  density  measured  by  the  Sondrestrom 
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Figure  4.  Simulated  radar  elevation  scans  through  the  simulated  volume  of  section  3.1.  (a)  The  large-scale  density  structure 
corresponds  to  the  TOl.  (b)  -  (i)  The  fragmentation  of  the  TOI  and  formation  of  smaller-scale  structures. 
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radar  -12  min  after  the  start  of  the  break-off  event.  Each 
figure  displays  two  density  structures  separated  by  a 
channel  of  low  density  located  almost  overhead.  The 
peak  density  in  each  of  the  structures  shows  a  good 
agreement  between  modeled  and  measured  values. 
However,  the  minimum  density  situated  between 
structures  is  much  more  depleted  in  the  experimental 
data.  It  is  below  105  cm'3  in  the  radar  data  but  above 
2xl05  cm'3  in  the  model.  Other  differences  between  the 
measured  and  modeled  densities  are  (1)  the  larger  north- 
south  extension  of  the  density  depleted  region  of  the 
measured  data,  and  (2)  the  noticeable  higher  bottomside 
F  layer  density  of  the  experimental  data. 

While  the  stationary  vortices  scheme  increases  the 
amount  of  density  structures  transiting  within  the  polar 
cap,  it  fails  to  produce  isolated  islands  of  enhanced 
density  structures  detached  from  the  auroral  oval.  A 
camera  located  within  the  polar  cap  will  be  able  to 
indicate  this  effect.  Moreover,  our  definition  of  polar  cap 
patches  postulates  that  less  dense  plasma  should  surround 
the  higher  density  plasma  in  all  directions  in  order  to 
qualify  as  a  patch.  Thus  we  conclude  that  stationary 
vortices  do  not  reproduce  all  the  characteristics  of  polar 
cap  patches. 

3.2.  Traveling  Vortices 

The  second  modeling  effort  was  carried  out  using  a  set 
of  two  traveling  vortices  moving  with  the  background 
plasma  flow.  To  generate  the  motion  of  the  vortices,  the 
two-vortex  system  was  displaced  in  the  antisunward 
direction  as  a  whole  entity.  The  displacement  was  set  by 
the  velocity  of  the  plasma  at  the  center  of  the  line  joining 
the  vortices’  centers.  Similar  to  the  previous  case,  the 
vortices  potential  was  introduced  at  1300  UT.  The 
potential  was  then  reduced  by  half  15  min  later.  The 
location  of  the  vortices  in  the  polar  cap  was  updated 
every  30  s;  this  interval  was  equal  to  the  step  time  of  the 
simulation.  Figure  5a  repeats  the  initial  density  of  Plate  1. 
During  the  first  9  min  of  the  simulation  (Figures  5b  -  5d), 
the  TOI  is  drastically  distorted  in  a  fashion  similar  to  the 
nonconvective  case.  However,  the  width  of  the  depletion 
channel  and  the  depth  of  the  density  decrease  are  larger 
for  the  convecting  vortices  case.  The  dawn-dusk 
extension  of  the  density  structure  is  smaller  in  this 
convective  scenario.  We  call  this  structure  a  region  of 
enhanced  density,  because  Ne  is  higher  inside  the 
structure  than  in  the  rest  of  the  polar  cap.  Figures  5d 
through  5f  show  that  both  density  structures  develop  an 
east-west  or  dawn-dusk  alignment.  They  are  also  well 
collocated  with  the  corresponding  vortices.  In  contrast  to 


the  nonconvecting  vortices  case,  we  see  (Figure  5f)  that 
the  more  northern  density  structure  has  no  attachment  to 
the  oval  or  subauroral  plasma.  The  last  three  panels  at  the 
bottom  of  Figure  5  show  that  as  two  patches  travel 
antisunward  across  the  polar  cap,  they  suffer  little 
distortion.  This  feature  is  also  different  from  the  previous 
case,  where  only  one  structure  was  fully  formed.  The 
prominent  characteristic  continues  to  be  that  both 
structures  are  individual  entities  with  no  attachment  to  the 
oval  or  subauroral  plasma.  The  polar  cap  patches,  as 
described  here,  are  very  similar  in  shape  to  the  “cigar” 
patches,  elongated  in  the  east-west  direction,  described  by 
Fukuietai  [1994]. 

A  more  quantitative  view  of  the  formation  of  the 
density  enhancement  and  depletion  is  displayed  in  Figure 
6.  Figures  6b  and  6c  show  the  cross  section  in  the 
magnetic  meridian  plane  of  two  structures  that  were 
formed  by  the  convective  vortices.  The  structure  located 
at  100  km  south  of  the  station  was  formed  from  plasma 
located  at  subauroral  latitudes  (south  of  the  TOI).  The  Ne 
structure  located  at  300  km  north  is  formed  from  the  TOI. 
The  particular  location  of  the  vortex  causes  plasma  from 
the  fragmented  TOI  to  intrude  into  the  southern  velocity 
structure  and  reach  values  near  9xl05  cm"3.  The  minimum 
density  between  patches  is  105  cm'3,  as  seen  at  100  km 
north  in  Figure  6e.  This  depletion  is  replenished  partially 
by  solar  radiation  and  reaches  values  near  3xl05  cm'3  in 
Figure  6h.  Comparison  of  Figures  6e  and  2c  indicates  a 
good  qualitative  agreement.  The  modeled  densities  of 
Figure  6e  show  a  slightly  higher  density  value  and  a 
narrower  width  in  the  region  of  depleted  densities.  In  this 
region,  the  modeled  density  are  <105  cm'3  at  300  km 
altitude,  similar  to  the  experimental  values. 

The  traveling  vortices  scenario  brought  only  a 
qualitative  agreement  with  the  density  of  Figure  2c. 
However,  it  was  able  to  produce  islands  of  high  density 
transiting  the  polar  cap.  Thus  we  conclude  that  the 
convecting  vortices  are  able  to  form  patches  of  enhanced 
densities  when  they  are  located  near  the  boundary 
between  the  polar  cap  and  the  auroral  oval. 

3.3.  Heppner-Maynard-Type  Bv  Negative  Globa! 
Pattern 

As  mentioned  in  section  2.1,  the  radar  data  is  not 
effective  at  discriminating  between  different  global 
patterns.  Thus  it  is  of  interest  to  see  how  our  results 
might  change  if  a  different  global  pattern  is  used  with  the 
vortices.  For  this  purpose,  we  repeated  the  calculations  of 
the  previous  section  to  determine  the  effect  on  patch 
formation  using  a  Heppner- Maynard-type  pattern.  Figure 
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Figure  5.  Same  as  Figure  3  but  for  the  modeling  study  of  traveling  convective  vortices  (section  3.2). 


7a  displays  the  initial  conditions  of  the  NmF2  values 
before  the  start  of  the  simulation.  The  TOI  shows  a  few 
peculiar  differences  with  respect  to  the  Heelis-type 
convection  pattern  of  section  3.1.  The  maximum  density 
is  slightly  higher,  and  the  TOI  is  a  much  wider  before  it 
turns  poleward.  These  differences  are  attributed  to  a 
slower  flow  velocity  (stagnation  region)  near  the  dayside 
region  in  the  Heppner-Maynard  global  convection 
patterns.  The  initial  location  of  the  vortices  and  the 
amplitude  of  each  vortex  potential  are  similar  to  the  case 
described  in  section  3.2. 


Figures  7b  and  7c  show  that  the  vortices  are  able  to 
bisect  a  region  of  the  TOI  and  orient  the  major  axis  of  the 
density  structure  in  the  dawn-dusk  direction.  A  second 
region  of  depleted  densities  is  seen  at  the  lower  latitudes 
in  Figures  7c  through  7f;  however,  this  structure  is  not 
completely  isolated  from  the  TOI.  Moreover,  it  joins  the 
new  TOI  that  is  being  formed  during  the  last  stages  of  the 
simulation.  The  patch  located  at  higher  latitudes  retains  its 
alignment  for  several  minutes  after  the  potential  was 
reduced  by  50%  (Figures  7g  -  7h).  Figure  8  shows  that 
the  patch  is  bounded  on  the  equatorward  side  by  plasma 
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of  traveling  convective  vortices  (section  3.2). 
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Figure  7.  Same  as  Figure  3  but  for  the  Heppner-Maynard  convection  pattern  (section  3.3). 


as  low  as  10s  cm'3.  It  also  demonstrates  that  the  patch 
moves  toward  the  pole  and  subsequently  breaks  from  the 
TOI.  When  the  patch  is  circulating  across  the  polar  cap, 
the  peak  density  inside  the  patch  is  8xl05  cm'3,  a  value 
much  higher  than  the  typical  l-2xl05  cm'3  seen  in  the 
polar  cap  and  outside  the  TOI.  Figures  8e  and  8f  show  a 
more  pronounced  depletion  between  the  density 
structures  as  compared  to  the  Heelis/traveling  vortices 
scenario.  When  compared  to  the  experimental  values  of 
Figure  2c,  we  notice  that  several  features  are  well 


reproduced.  The  altitude  extension  of  both  structures  are 
very  similar.  N,  =  7xl05  cm'3  at  700  km  is  observed  in 
both  experimental  and  modeled  data.  The  region  of 
depleted  density  is  wider  than  in  the  previous  case.  The 
poleward  motion  of  the  structures  also  resembles  more 
closely  the  movement  of  the  structures  seen  in  the 
experimental  data.  The  Heppner-Maynard  pattern 
successfully  produced  a  density  structure  with  all  the 
characteristics  of  a  polar  cap  patch.  In  contrast  to  the 
Heelis  convection  pattern,  it  creates  a  single  structure. 


78 


19  FEB  90.  .  .Density*. IQ'1  .  .  .1300UT  1 9. FEB  90  .  .  Density*! O'5  .  .  .1303  UT  19. FEB  90  .  Density*  10.5  ,  .  .1306  UT 


convection 


588 


Geomagnetic  Coordinates 


Figure  9.  Same  as  Figure  3  but  after  eliminating  the  dependence  of  the  ion  temperature  on  the  magnetospheric 
electric  field  (section  3.4). 


3.4.  Assessment  of  the  Temperature  Effect  on  Patch 
Formation 

The  fourth  case  study  was  performed  to  quantify  the 
importance  of  the  0+  recombination  rate  in  producing 
channels  of  low  density.  The  results  of  this  study  are 
presented  in  Figure  9.  In  this  series  of  simulations,  we 
returned  to  the  Heelis/traveling  vortices  scenario  except 
we  turned  off  the  dependence  of  Tt  on  the  magnetospheric 
electric  field  and  hence  turned  off  any  enhanced  0+  loss 


rates.  The  goal  in  this  investigation  was  to  determine 
quantitatively  how  much  plasma  depletion  is  attributed  to 
erosion  by  enhancements  of  the  0+  loss  rate  and  what 
percentage  of  the  plasma  depletion  is  carried  from  earlier 
local  times  by  the  vortex  flow.  Evidently,  an  eastward 
flow  is  needed  in  order  to  transport  the  less  dense 
morning  side  plasma  into  the  afternoon  region.  A  flow  of 
this  type  exists  at  the  common  region  between  the 
vortices. 
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dependence  on  the  magnetospheric  electric  field. 


Figure  11.  Same  as  Figure  3,  but  for  nine  different  initial  locations  of  the  traveling  vortices. 


A  comparison  of  Figures  9  and  5  indicates  some 
resemblance  between  both  figures.  The  structuring  is  also 
initiated  with  a  rotation  of  the  TOI  toward  a  dawn-dusk 
alignment.  Plasma  from  the  early  morning  side  also 
intrudes  equatorward  of  the  northern  patch.  However, 
there  are  also  some  prominent  differences.  We  did  not 
find  a  full  separation  between  the  Nc  structure  and  the 
auroral  oval,  even  15  min  into  the  simulation.  There  is  a 
channel  containing  plasma  density  equal  to  4xl05  cm'3  in 
the  afternoon  side  of  the  structure.  The  southern  Ne 


structure  is  absent,  because  at  this  latitude  the  southern 
vortex  can  only  bring  more  dense  plasma  from  the 
afternoon  sector.  Figure  I  Of  indicates  that  the  deepest 
depletion  achieved  in  this  simulation  is  4x10s  cm'3,  while 
the  minimum  density  value  during  the  full-blown 
simulation  was  10s  cm 3  in  Figure  5.  Considering  that  the 
original  density  was  9xl05  cm'3,  then  62%  of  the 
depletion  is  due  to  convection  of  low-density  plasma  and 
38%  is  due  to  erosion  by  the  temperature-dependent  loss 
rate. 
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4.  Traveling  Vortices  at  Different  Local 
Times 

In  the  previous  section,  we  were  concerned  with  the 
formation  of  patches  under  different  patterns  of  global 
convection.  In  this  section,  we  search  for  the  preferred 
local  time  at  which  the  vortices  are  more  likely  to 
generate  polar  cap  patches. 

Traveling  vortices  have  been  observed  at  a  large 
variety  of  different  local  times  [Friis-Christensen  et  al ., 
1988;  Heikkila  et  al ,  1989]  using  radar  and 

magnetometer  data.  They  have  also  been  predicted  on 
theoretical  grounds  [Newell  and  Sibeck ,  1993].  They 
postulated  an  association  between  the  presence  of  twin 
vortices  and  changes  in  the  high-latitude  convection 
patterns  due  to  pressure  pulses  in  the  magnetosphere.  In 
order  for  the  vortices  to  produce  patches,  evidently,  they 
have  to  develop  near  the  TOI  and  have  a  lifetime  long 
enough  to  erode  a  sector  of  the  TOI.  We  have  carried  out 
eight  additional  simulations  placing  the  vortices  at 
different  local  times.  The  local  times  that  were  used 
varied  from  0800  to  1600  LT.  The  results  of  this  study  are 
compiled  in  Figure  1 1.  This  figure  compares  NmF2  values 
obtained  15  min  along  the  simulations  (at  1315  UT)  for 
the  eight  new  cases  and  the  simulation  described  in 
section  3.2.  The  number  on  the  upper  left  comer  indicates 
the  local  times  at  which  the  vortices  were  located  at  the 
beginning  of  the  simulation.  Only  the  vortices  located  at 
1000  LT,  1100  LT  and  1200  LT  (Figures  11c,  lid  and 
lie)  produce  distinct  patches.  Vortices  at  0800  LT  and 
0900  LT  (Figures  1  la  and  1  lb)  only  distort  the  TOI  and 
create  isolated  regions  of  low  density  (2xl05  cm'3) 
plasma.  When  these  structures  circulate  within  the  polar 
cap,  they  will  consist  of  plasma  equally  dense  to  the  polar 
cap  plasma.  When  the  vortices  were  located  at  1300  LT 
(Figure  Ilf),  we  obtained  dawn-dusk  elongated  density 
traces.  No  detachment  from  the  auroral  oval  densities 
was  observed.  The  last  three  modeling  efforts  in  this 
subsection  were  performed  using  vortices  located  at  1400 
LT,  1500  LT  and  1600  LT  (Figures  llg-lli).  These 
figures  depict  the  formation  of  some  kind  of  structuring  in 
the  afternoon  cell.  These  structures  are  connected  to  the 
auroral  oval  and  will  not  penetrate  into  the  polar  cap.  In 
fact,  we  followed  all  these  case  simulations  for  another  30 
min  and  found  the  structures  continuously  elongating  in 
the  dawn-dusk  direction.  These  structures  could  certainly 
appear  as  isolated  structures  to  digisondes  or  to  low- 
orbiting  satellites,  but  an  all-sky  camera  properly  located 
could  possibly  indicate  the  connection  of  these  structures 
to  the  auroral  plasma. 


5.  Discussion  and  Conclusions 

Several  researchers  have  indicated  that  the  TOI  is 
present  at  longitudes  near  Sondrestrom  when  this  high- 
latitude  station  is  near  noon  [Foster  and  Doupnik,  1984; 
Kelly  and  Vickrey ,  1984].  Our  simulation  during 
undisturbed  conditions  has  also  demonstrated  that  the 
TOI  is  a  prominent  feature  at  Sondrestrom  when  we 
select  a  typical  Bz  south  convection  pattern.  The 
mechanism  described  here  requires  the  existence  of  an 
elongated  TOI  as  the  basic  ingredient  to  form  a  patch.  It 
also  needs  a  large  plasma  jet  occurring  near  and  across 
the  TOI.  The  large  plasma  velocity  will  create  a  region  of 
enhanced  Th  consequently  increase  the  0+  recombination 
rate  and  then  erode  the  local  plasma.  We  found  also  that 
low-density  plasma  can  be  brought  in  by  the  vortex  flow 
from  earlier  local  times.  This  patch  formation  mechanism 
will  create,  in  this  way,  sections  of  depleted  and  enhanced 
densities  across  the  TOI. 

The  elevation  scans  conducted  across  our  simulated 
volume  demonstrated  that  the  TOI  or  other  structures  in 
the  polar  cap  could  be  wrongly  interpreted  as  isolated 
entities  (say  patches).  In  reality,  they  may  be  connected 
to  the  auroral/subaurora!  plasma.  We  conclude  that  it  is 
necessary  to  conduct  antenna  azimuthal  scans  or  to  use  an 
all-sky  camera  to  properly  identify  a  polar  cap  patch. 

Sojka  et  al  [1993,  1994]  have  used  a  scheme  of 
continuous  changes  of  the  By  component  of  the  IMF  to 
redirect  the  TOI  toward  the  dawn  or  dusk  sides  of  the 
polar  cap.  In  this  study,  we  have  employed  the 
appearance  of  vortices  in  the  polar  cap  and  their 
subsequent  motion  across  the  polar  cap  to  produce  polar 
cap  patches.  No  repetitive  changes  in  any  of  the 
components  of  the  IMF  were  required.  The  proposed 
mechanism  only  needs  the  existence  of  vortices  moving 
with  the  background  plasma.  Moreover,  it  was  found  that 
the  convecting  vortices  scenario  was  the  most  favorable 
situation  to  form  disconnected  plasma  structures  inside 
the  polar  cap. 

The  modeled  densities  of  Figures  6  and  8  indicated  a 
good  qualitative  agreement  with  the  data  measured  by  the 
Sondrestrom  ISR.  The  discrepancies  between  measured 
and  modeled  densities  can  be  understood  in  terms  of  the 
initial  inputs  that  were  used  in  the  simulations  and  a  few 
of  the  limitations  of  the  GTIM  model.  The  much  smaller 
density,  seen  in  the  depleted  density  region  of  Figure  2c, 
can  be  attributed  to  a  much  smaller  density  in  the 
morning  cell  from  where  the  vortices  grab  the  low- 
density  plasma.  In  fact,  the  Sondrestrom  ISR  detected 
values  near  2xl05  cm'3  at  1100  UT  (2  hours  before  the 
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measurements  reported  here)  that  if  unchanged  could 
explain  this  discrepancy.  The  radar  also  measured 
densities  above  105  cm4  at  200  km  altitude;  no  similar 
feature  was  reported  in  the  simulations  because  of  the 
lack  of  molecular  ions  in  our  simulations.  The  GTIM 


model  is  at  the  present  time  essentially  a  one-ion  (0+) 
model.  Simulations  were  performed  with  and  without 
particle  precipitation  in  the  auroral  oval,  and  essentially 
no  differences  were  seen  in  the  F  region.  The  effects  of 
particle  precipitation  associated  with  the  vortex  located 
more  equatorward  where  convergent  electric  fields 
existed  were  not  included  and  neither  were  the  effects  of 
soft  precipitation  that  can  occur  in  the  cusp/cleft  and  polar 
cap.  These  issues  will  be  topics  for  future  work. 

We  found  that  62%  of  the  density  depletion  was  due  to 
plasma  being  carried  from  earlier  local  times.  This  fact 
argues  in  favor  of  the  ability  of  smaller  potentials 
associated  with  the  vortices  to  also  create  regions  of 
depleted  density.  However,  the  vortex  potential  cannot 
be  too  small,  otherwise  the  vortex  velocity  will  also  be 
small,  and  the  time  to  transport  low-density  plasma  from 
the  morning  cell  will  be  considerably  longer. 

From  this  study  we  found  the  following: 

1.  Polar  cap  patches  can  be  generated  by  traveling 
vortices  independent  of  the  convection  pattern  being 
used.  We  have  used  a  Heelis-type  and  a  Heppner- 
Maynard-type  convection  pattern  for  By  positive  and  for 
By  negative,  respectively.  An  essential  condition  for 
forming  well-separated  patches  was  to  allow  the  vortices 
to  travel  with  the  background  global  convection. 

2.  Enhanced  recombination  of  0+  contributes  to  the 
creation  of  a  channel  of  low  density.  Equally  important  is 
the  transport  of  low-density  plasma  from  earlier  local 
times.  For  the  size  and  potential  of  the  vortices  that  we 
used  38%  of  the  depletion  was  attributed  to  the  0+ 
recombination  loss  and  62%  to  transport. 

3.  There  is  a  preferential  local  time  at  which  the 
vortices  can  generate  patches.  We  found  that  this  local 
time  sector  is  restricted  to  between  1000  LT  and  1200  LT. 


4.  The  simulation  presented  here  qualitatively  agrees 
with  the  data  collected  at  Sondrestrom  on  February  19, 
1990  [Valladares  et  al. ,  1994b];  both  exhibit  the  same 
salient  features. 

5.  Our  modeling  described  here  postulates  that  the 
polar  cap  patches  will  intrinsically  have  the  shape  of  the 
vortices.  Circular  vortices  will  reproduce  more  circular 
patches.  Elliptical  vortices,  as  presented  here,  will 
generate  the  elongated  cigar-shaped  patches  that  have 
been  found  in  the  polar  cap  [Fukui  et  al. ,  1 994]. 
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Collisional  degradation  of  the  proton-H  atom  fluxes  in 
the  atmosphere:  A  comparison  of  theoretical 
techniques 

Dwight  T.  Decker,1  Boris  V.  Kozelov,2  B.  Basu,3  J.  R.  Jasperse,3 
and  V.  E.  Ivanov2 


Abstract.  Three  methods  for  calculating  the  transport  of  energetic  protons 
and  hydrogen  atoms  within  the  Earth’s  atmosphere  are  compared.  The  methods 
are  (1)  a  Monte  Carlo  (MC)  simulation,  (2)  a  discrete  energy  loss  solution  to 
the  linear  transport  equations,  and  (3)  a  continuous  slowing-down  approximation 
(CSDA).  In  the  calculations  performed,  all  three  models  use  the  same  cross  sections, 
three-component  (N2,02,0)  neutral  atmosphere,  and  incident  isotropic  Maxwellian 
proton  fluxes  of  various  characteristic  energies  (1-20  keV).  To  ensure  that  all  three 
methods  include  the  same  physical  processess,  the  effects  of  magnetic  mirroring 
and  the  lateral  spreading  of  particles  are  “turned  off”  in  the  MC  simulations  as 
these  processes  are  not  included  in  the  present  linear  transport  or  CSDA  models. 
A  variety  of  quantities  are  calculated  and  compared  including  energy  deposition 
rates,  eV/ion  pair,  hemispherically  averaged  differential  fluxes  of  protons  and  H 
atoms,  energy  integrated  differential  fluxes,  and  total  proton  and  H  atom  fluxes. 
The  agreement  between  all  three  models  is  excellent  except  at  the  lowest  altitudes. 
Apart  from  these  altitudes,  the  differences  that  do  exist  are  small  compared  to  the 
errors  that  generally  result  from  poorly  known  inputs  and  compared  to  the  typical 
errors  quoted  for  geophysical  observations.  The  altitudes  where  the  results  do  differ 
significantly  are  where  the  proton  and  H  atom  fluxes  are  severely  attenuated  and 
are  below  the  altitudes  where  the  bulk  of  the  energy  deposition  and  ionization  takes 
place.  The  success  of  these  comparisons  suggests  that  our  ability  to  model  actual 
observations  is  presently  limited  by  uncertainties  in  cross  sections  and  the  lack  of 
suitable  observations  rather  than  our  ability  to  solve  the  equations  that  describe 
the  known  physics  of  proton-H  atom  transport. 


1.  Introduction 

It  is  well  established  that  an  essentially  permanent 
feature  of  the  Earth’s  high-latitude  atmosphere  is  the 
presence  of  auroral  particle  fluxes.  It  is  also  well  known 
that  these  precipitating  fluxes  from  the  magnetosphere 
consist  of  mostly  electrons  and  protons  with  a  small 
admixture  of  other  ions.  Early  ground-based  optical 
observations  [Romick  and  Elvey,  1958;  Galperin,  1959] 
indicated  that  electron  and  proton  aurora  could  have 
separate  latitudinal  distributions.  Later  satellite  obser¬ 
vations  of  the  precipitating  particles  established  that 
statistically  the  electron  and  protons  precipitate  within 
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two  ovals  that  are  not  colocated  [Shorter,  1981;  Hardy 
et  al.,  1989].  Further,  an  examination  of  the  Hardy  sta¬ 
tistical  models  [Hardy  et  al,  1987,  1991]  reveals  that  in 
the  afternoon  and  evening  sectors  the  ions  at  the  lower 
auroral  latitudes  can  carry  a  significant  portion  of  the 
incoming  energy  flux.  Thus,  in  order  to  study  the  dissi¬ 
pation  of  auroral  energy  within  the  Earth’s  atmosphere, 
it  is  necessary  to  study  not  just  the  transport  and  colli¬ 
sional  degradation  of  energetic  electrons  but  also  those 
of  protons. 

The  transport  of  energetic  protons  within  matter  is  a 
problem  that  has  been  studied  for  a  variety  of  physical 
situations.  In  the  case  of  transport  of  auroral  protons 
within  the  atmosphere  there  are  several  features  that 
add  to  the  difficulty  of  the  problem.  Some  of  those  fea¬ 
tures  are  (1)  auroral  proton  energies  can  be  less  than 
the  energy  of  the  maximum  of  the  ionization  cross  sec¬ 
tion,  (2)  there  is  a  large  gradient  of  the  atmospheric 
density  with  respect  to  altitude,  and  (3)  the  Earth’s 
magnetic  field  is  present.  The  first  feature  requires 
detailed  knowledge  of  cross  sections  since  simple  esti¬ 
mations  (Born  approximation)  of  cross  sections  are  not 
/alid.  Further,  at  auroral  energies  the  charge  chang- 
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ing  processes  have  to  be  included  and  this  introduces 
energetic  hydrogen  (H)  atoms  into  the  problem.  The 
effect  of  the  large  gradient  in  the  atmosphere  is  to  cre¬ 
ate  a  strong  altitude  dependence  in  the  mean  free  path 
of  the  protons  and  H  atoms.  Hence  the  behavior  of 
the  particle  flux  is  highly  dependent  on  altitude.  The 
presence  of  the  Earth’s  magnetic  field  acts  to  cause  the 
protons  to  follow  field  lines  into  the  ionosphere.  How¬ 
ever,  the  H  atoms  are  not  constrained  by  the  magnetic 
field.  This  coupled  with  the  altitude-dependent  mean- 
free  path  results  in  a  lateral  spreading  of  the  p-H  flux 
at  high  altitudes  but  essentially  no  spreading  at  low  al¬ 
titudes.  Finally,  a  feature  that  can  be  a  useful  simplifi¬ 
cation  in  any  modeling  is  that  all  the  collision  processes 
involving  protons  and  H  atoms  are  strongly  peaked  in 
the  forward  direction.  However,  at  very  low  energies 
(E  <  1  keV)  energy  loss  due  to  elastic  collisions  does 
become  significant. 

There  have  been  several  treatments  of  p-H  transport 
that  have  appeared  in  the  literature.  They  can  be  gen¬ 
erally  categorized  as  range  theoretic,  continuous  slow¬ 
ing  down,  linear  transport  theoretic,  and  Monte  Carlo. 
The  range  theoretic  method  is  based  on  laboratory  mea¬ 
surements  of  range  energy  relations  for  protons  in  air 
or  its  theoretical  approximations  and  has  been  used  in 
the  works  of  Eather  and  Burrows  [1966],  Eather  [1967, 
1970],  Isaev  and  Pudovkin  [1972],  Eenriksen  [1979],  and 
Rees  [1982].  The  continuous  slowing-down  approxima¬ 
tion  (CSDA)  uses  more  detailed  information  about  scat¬ 
tering  processes  and  was  used  for  the  auroral  proton 
problem  by  Edgar  et  al.  [1973,  1975],  Work  on  analytic 
estimations  of  lateral  spreading  of  p-H  flux  [ Johnstone , 
1972;  Iglesias  and  Vondrak ,  1974]  also  can  be  classified 
as  using  the  CSDA. 

In  more  recent  years,  linear  transport  (LT)  theory  has 
been  used  to  solve  for  the  auroral  proton  and  H  atom 
fluxes.  In  the  work  of  J asperse  and  Basu  [1982]  the 
coupled  transport  equations  based  on  the  Boltzmann 
equation  were  solved  analytically,  and  in  a  follow-up 
paper  [Basu  et  al .,  1987],  theoretically  calculated  elec¬ 
tron  density  profiles  due  to  incident  protons  were  found 
in  good  agreement  with  Chatanika  radar  data.  This  in 
turn  was  followed  by  the  development  of  a  fully  numeri¬ 
cal  code  for  solving  the  linear  transport  equations  [Basu 
et  a/.,  1990]. 

The  Monte  Carlo  (MC)  method  has  been  applied  to 
the  auroral  proton  problem  by  Davidson  [1965],  Pono¬ 
marev  [1976],  Galperin  et  al.  [1976],  and  Kozelov  and 
Ivanov  [1992].  Most  recently,  Kozelov  [1993]  has  imple¬ 
mented  the  method  using  a  full  three  species  neutral 
atmosphere  and  the  influence  of  the  Earth’s  magnetic 
field.  While  computationally  intensive  the  Monte  Carlo 
method  does  have  the  advantage  of  including  processes, 
such  as  beam  spreading  or  mirroring,  that  are  difficult 
to  treat  in  the  other  approaches. 

Assessing  the  validity  of  any  approach  usually  in¬ 
volves  either  “back  of  the  envelope”  estimates  of  ne¬ 
glected  processes,  a  posteriori  calculations  of  those  pro¬ 
cesses,  comparisons  to  observations,  or  comparisons  to 
alternative  approaches.  The  difficulty  in  comparing  ap¬ 
proaches  is  that  when  comparing  to  results  of  other 


authors  there  are  often  differences  in  cross  sections, 
boundary  conditions,  and  geophysical  conditions  that 
confuse  the  issue  of  comparing  methods.  MC  and  CSDA 
models  with  the  same  input  parameters  were  compared 
by  Porter  and  Green  [1975],  but  the  comparison  was 
only  for  homogeneous  molecular  nitrogen,  and  the  ini¬ 
tial  proton  flux  was  taken  as  isotropic  and  monoener- 
getic  with  initial  energy  1  keV.  More  extensive  compar¬ 
isons  for  a  N2  atmosphere  were  made  by  Kozelov  and 
Ivanov  [1992].  They  considered  initial  energies  form  1 
keV  to  1  MeV.  For  initial  energies  less  than  10  keV, 
differences  were  found  near  the  high-altitude  boundary 
and  at  lowest  altitudes  of  penetration. 

The  purpose  of  this  paper  is  to  compare  three  meth¬ 
ods  for  calculating  the  transport  and  coilisional  degra¬ 
dation  of  energetic  protons  and  H  atoms  in  the  Earth’s 
ionosphere  using  identical  cross  sections,  boundary  con¬ 
ditions,  and  geophysical  conditions.  We  will  show  that 
for  several  quantities  the  differences  between  the  three 
methods  are  reasonably  small  compared  to  the  model 
uncertainties  that  arise  from  uncertainties  in  the  re¬ 
quired  inputs  (cross  sections,  neutral  densities,  and 
boundary  conditions)  and  the  uncertainties  in  the  avail¬ 
able  observations  (particle  fluxes,  optical  emissions,  and 
ion  densities).  The  three  methods  considered  are  the 
CSDA,  LT,  and  MC.  In  section  2  we  briefly  describe 
the  three  models  used.  A  more  detailed  discussion  of 
the  relationship  between  the  CSDA  and  LT  models  is 
given  in  an  appendix.  After  a  discussion  in  section  3 
concerning  the  physical  processes  and  input  parame¬ 
ters  included  in  all  three  models,  we  present  in  section 
4  the  comparison  of  the  various  results  from  the  models. 
The  paper  concludes  with  a  discussion  and  summary  in 
section  5. 


2.  Description  of  Theoretical  Models 


2.1.  Continuous  Slowing-Down  Approximation 

The  continuous  slowing-down  approximation  (CSDA) 
describes  the  degradation  of  the  energetic  particle’s  en¬ 
ergy  in  terms  of  loss  functions.  For  a  monoenergetic 
monodirectional  p-H  flux  the  equation  for  energy  loss  is 
given  by 


dE  /  \ 


La,p(J5)$p(z,£,^)  ^ 


,-E,/*)] 


$p(z,E,fi)  +  $B(z,E,(i) 


(1) 


where  $  values  are  the  energetic  particle  fluxes,  z  is  the 
altitude,  p  is  cosine  of  the  pitch  angle,  and  na(z)  are 
the  densities  for  neutral  specie  a.  LatP(E)  and  La^js(E) 
are  the  loss  functions  for  the  protons  and  H  atoms  due 
to  collisions  with  neutral  gases  and  are  given  by 

L*AE)  =  £<p(£)<p(£)  (2) 

3 

'  La<B(E)  =  £<*(£)<*(£)  (3) 
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where  cr^  P{E)  is  the  cross  section  for  a  type  j  collision 
between  a  proton  and  a  neutral  specie  a  and  Wl  p(E)  is 
the  corresponding  energy  loss.  Similarly,  o-’a  ff(E)  and 

are  the  cross  sections  and  the  energy  loss  for 
collisions  involving  H  atoms. 

Usually,  the  loss  function  is  defined  for  a  charge  state 
equilibrium  flux  [Porter  and  Green ,  1975].  Here,  in  an 
implementation  developed  by  one  of  us  (B.V.K.),  the 
loss  function  is  defined  using  the  nonequilibrium  flux 
that  describes  the  modification  of  the  flux  charge  state 
composition.  The  nonequilibrium  flux  is  calculated  us¬ 
ing  the  following  equations: 

a 

-  $P{z, U, n)  na(2)<ri°(J5)  (4) 

_  d$P(z,E,  v) 
dz  dz 

where  <r10  and  a01  are  the  total  cross  sections  for  charge 
exchange  and  stripping  respectively.  The  calculational 
procedure  consists  of  solving  (1),  (4),  and  (5)  on  a  alti¬ 
tude,  energy  and  pitch  angle  grid  with  the  incoming  flux 
at  the  high-altitude  boundary  being  approximated  by  a 
set  of  monoenergetic  monodirectional  particle  streams. 
A  more  detailed  discussion  of  the  CSDA  can  be  found 
in  the  appendix. 

2.2.  Linear  Transport  Model 

The  fundamental  feature  of  the  linear  transport  (LT) 
model  is  the  use  of  transport  theoretic  methods  to  solve 
two  coupled  Boltzmann  equations  for  the  proton  and 
H  atom  differential  fluxes.  In  solving  these  equations, 
we  assume  steady  state  conditions,  no  electric  fields, 
a  uniform  geomagnetic  field  (no  magnetic  mirroring), 
and  azimuthal  symmetry  about  the  geomagnetic  field 
line.  The  collisional  processes  included  are  the  various 
types  of  collisions  between  the  energetic  precipitating 
particles  and  the  neutral  constituents  of  the  ionosphere. 
The  protons  and  H  atoms  are  coupled  to  each  other 
by  charge-changing  collisions  with  neutrals  (charge  ex¬ 
change  and  stripping).  The  two  transport  equations, 
discussions  of  the  various  collision  processes,  and  the 
numerical  methods  for  solving  the  equations  can  be 
found  in  the  work  of  Basu  et  tal  [1993]  and  will  not 
be  given  here.  A  discussion  of  the  relationship  between 
the  linear  transport  theory  and  the  CSDA  can  be  found 
in  the  appendix. 

2.3.  Monte  Carlo  Model 

The  Monte  Carlo  (MC)  model  uses  a  “collision-by¬ 
collision”  algorithm  which  is  based  on  simulating  the 
individual  trajectories  of  a  large  number  of  protons  and 
H  atoms.  The  “history”  of  each  particle  is  represented 
as  a  series  of  collisions  that  are  separated  by  segments 
of  free  streaming  within  a  magnetic  field.  The  selection 
of  the  path  length  between  collisions  and  the  collisional 
parameters  is  made  by  mapping  into  a  random  number 
f  generated  uniformly  in  the  range  0-1. 


The  altitude  of  the  particle  after  the  (t  -f  l)5*  path 
length  is  determined  from  the  formula 

*<+i 

-ln*=  f  7E"»WffS(^  (6) 
i  a 

where  z{  is  the  altitude  before  the  (i  +  l)s‘  path  length, 
caj}(Ei)  is  the  total  cross  section  of  the  particle  in  gas 
q,  E{  is  the  particle’s  energy  after  the  »th  collision.  /3 
is  the  charge  state  of  the  particle  (p  or  H).  The  neu¬ 
tral  species  with  which  the  ( i  -f  l)5t  collision  occurred 
are  determined  by  random  sampling  using  the  proba¬ 
bilities  ntt(z)/X^na(z),  where  a  is  N2,  02)  or  O.  The 
type  of  collision  is  also  selected  by  random  sampling 
using  probabilities  a’(Ei)/atoi(Ei),  where  j  is  the  type 
of  collision.  When  a  collision  occurs,  the  energy  of  the 
particle  is  decreased  based  on  the  type  of  collision  used. 
In  the  case  of  charge  exchange  or  stripping  collisions 
the  charge  state  of  the  particle  is  modified.  To  sim¬ 
ulate  an  initial  Maxwellian  energy  distribution  at  the 
high-altitude  boundary,  we  divide  the  full  energy  range 
into  10-12  energy  channels  and  place  within  each  chan¬ 
nel  8000  to  24,000  protons.  The  cosine  of  the  initial 
pitch  angle  for  each  particle  is  determined  from  the  ex¬ 
pression  /z  =  cos  0  =  \/£,  where  £  is  a  random  number 
generated  uniformly  in  the  range  0-1.  During  the  sim¬ 
ulation  of  a  particle’s  trajectory  its  location  is  stored 
along  with  the  various  quantities  needed  for  determin¬ 
ing  such  quantities  as  the  particle  fluxes,  energy  deposi¬ 
tion,  and  ionization  rates.  The  results  of  the  simulation 
are  processed  by  statistical  methods.  If  we  have  N  real¬ 
izations  of  a  random  variable  x,  then  the  average  value 

is  taken  as  x  =  (1/N)  x{.  The  statistical  error  Ax 

is  estimated  from  expressions 


Ax  =  C\/5x, 


where  x  is  the  random  value,  x*  is  the  ith  realization 
of  the  value,  Dx  is  the  dispersion.  If  C  &  1.4,  then  the 
probability  that  the  value  x  €  [x  —  Ax,  x  +  Ax]  is  equal 
to  0.7. 


3.  Inputs 

In  order  to  facilitate  our  comparison,  we  wanted  to 
make  the  geophysical  inputs,  the  microscopic  parame¬ 
ters,  and  the  physical  processes  included  in  the  three 
models  as  identical  as  possible.  Considering  the  physi¬ 
cal  processes,  we  find  that  all  three  models  include  the 
transport  and  degradation  of  protons  and  neutral  hy¬ 
drogen  atoms  as  they  collide  with  the  neutral  species  of 
the  ionosphere.  On  one  hand,  the  CSDA  and  LT  assume 
a  uniform  magnetic  field  (no  mirroring)  and  a  plane  par¬ 
allel  geometry  (no  spreading),  on  the  other  hand  the 
MC  method  can  handle  the  effects  of  magnetic  mirror¬ 
ing  as  well  as  lateral  spreading  of  the  energetic  particles. 
For  these  comparisons,  magnetic  mirroring  and  lateral 
spreading  were  turned  off  in  the  MC  calculations.  Sim¬ 
ilarly,  both  CSDA  and  LT  assume  forward  scattering, 
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while  MC  can  include  angular  scattering.  Again,  we 
simplify  the  MC  model  and  use  the  forward  scattering 
approximation  in  all  calculations  for  this  paper. 

By  microscopic  parameter  we  are  referring  to  the  var¬ 
ious  cross  sections  and  energy  losses  that  are  needed  by 
the  models.  The  MC  and  CSDA  use  a  more  extensive 
set  than  is  used  by  the  LT  model.  While  what  cross 
sections  are  used  can  have  serious  impact  when  compar¬ 
ing  to  observations,  the  critical  point  here  is  to  have  all 
three  models  use  the  same  set.  We  chose  to  use  the  sim¬ 
pler  LT  set  as  described  by  Basu  et  al.  [1993].  The  one 
exception  is  that  for  the  average  energy  of  the  ejected 
electron  in  the  stripping  process  we  use  the  same  ex¬ 
pression  as  is  used  for  the  ionization  process.  This  was 
done  simply  as  a  convenience  for  making  the  compari¬ 
son  between  the  three  models. 

Finally,  let  us  turn  to  what  we  call  the  geophysical 
conditions,  that  is,  the  neutral  atmosphere  and  the  in¬ 
coming  particle  flux.  The  neutral  atmospheric  densities 
(N2,  O2,  and  O)  and  the  neutral  temperature  as  a  func¬ 
tion  of  altitude  are  obtained  from  the  mass  spectrom¬ 
eter/incoherent  scatter  (MSIS-86)  neutral  atmosphere 
model  [ Hedin ,  1987].  The  MSIS-86  parameters  used 
were  an  F10.7  and  a  81-day  averaged  F10.7  value  of 
150,  a  dally  Ap  of  20,  a  geographic  latitude  of  65°  N,  a 
longitude  of  35°  E,  a  local  solar  time  of  24,  and  a  date 
of  December  16, 1993.  All  calculations  have  a  boundary 
altitude  at  700  km  and  an  incident  isotropic  proton  flux 
that  has  a  Maxwellian  energy  distribution. 

4.  Results 

Calculations  were  made  using  the  models  for  a  range 
of  Maxwellian  characteristic  energies  (1-20  keV)  and  an 
incident  energy  flux  of  0.5  ergs  cm"2s~1.  The  results  in 
all  cases  were  similar  so  we  will  present  detailed  com¬ 
parisons  for  just  one  case:  a  characteristic  energy  of  8 


keV  and  a  minimum  energy  of 500  eV.  In  Figures  la-lf 
,  we  show  the  hemispherically  averaged  differential  flux 
for  the  protons  and  H  atoms  at  several  selected  alti¬ 
tudes.  In  this  and  in  the  following  figures,  the  CSDA 
results  are  given  by  the  dotted  curves,  the  LT  by  the 
dashed  curves  and  the  MC  by  either  squares  or  error 
bars.  The  statistical  error  in  the  MC  results  is  either 
indicated  by  the  error  bars  or  is  smaller  than  the  size 
of  the  square  symbols  used  in  the  plots.  At  550  km, 
150  km  below  the  boundary  altitude,  we  see  that  the 
proton  energy  distributions  from  all  three  models  are  in 
excellent  agreement.  Likewise,  the  H  atom  plot  shows 
excellent  agreement,  though  the  MC  simulation  is  a  bit 
”  noisy”  because  of  the  lack  of  particles  as  the  H  atom 
flux  builds  up  from  a  starting  value  of  zero  at  the  bound¬ 
ary.  At  250  km  (Figure  lb)  the  agreement  remains  ex¬ 
cellent  between  all  three  models.  While  the  shape  of 
the  spectra  are  little  changed,  the  H  atom  flux  has  con¬ 
tinued  to  increase  at  the  expense  of  the  proton  flux. 
Descending  further  to  152  km  (Figure  lc),  we  see  that 
the  low-energy  part  of  the  spectra  have  filled  in  as  parti¬ 
cles  cascade  down  from  higher  energies.  The  MC  shows 
a  little  more  noise  at  the  lowest  energies,  but  the  agree¬ 
ment  remains  fairly  good  between  all  models.  Dropping 
to  lower  altitudes,  the  low-energy  portion  of  the  spec¬ 
tra  continue  to  increase  and  with  the  three  models  still 
in  reasonable  agreement  reach  a  maximum  around  118 
km  (Figure  Id).  Though  we  do  see  a  clear  separation 
between  the  CSDA  and  LT  results  at  lower  energies. 
At  110  km  (Figure  le),  the  magnitude  of  the  spectra 
at  all  energies  have  decreased  and  the  CSDA  and  LT 
continue  to  show  a  clear  separation.  Finally,  just  4  km 
lower  at  106  km  (Figure  If),  the  CSDA  and  LT  are 
showing  larger  differences,  and  the  error  in  the  MC  re¬ 
sults  have  grown  to  the  point  that  it  is  difficult  to  decide 
if  it  agrees  more  with  the  CSDA  or  the  LT  model. 

Having  considered  the  energy  distribution  of  the  par- 


Figure  la.  Hemispherically  averaged  (left)  proton  and  (right)  hydrogen  atom  fluxes  versus 
energy  at  an  altitude  of  550  km.  The  error  bars  are  the  Monte  Carlo  results,  the  long  dashes  are 
the  linear  transport  results,  and  the  dotted  curves  are  the  results  from  the  continuous  slowing- 
down  approximation.  The  incident  proton  flux  is  a  Maxwellian  with  a  characteristic  energy  of  8 
keV  and  a  total  incident  energy  flux  of  0.5  ergs  cm-^"1. 
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Figure  le.  Same  as  Figure  la  for  an  altitude  of  110  km. 


Figure  If,  Same  as  Figure  la  for  an  altitude  of  106  km. 


tide  fluxes  averaged  over  pitch  angle,  let  us  now  turn  to 
plots  (  Figures  2a-2g)  of  the  differential  flux  integrated 
over  energy  versus  the  cosine  of  the  pitch  angle.  At  670 
km  (Figure  2a),  we  see  an  isotropic  proton  flux  with  all 
three  models  in  excellent  agreement  and  a  much  smaller 
H  atom  flux  that  has  had  only  30  km  to  build  up  from 
zero  flux  at  the  boundary.  The  H  atom  flux  is  largest 
near  90°  pitch  angle  where  the  protons  have  undergone 
the  greatest  number  of  collisions.  Figures  2b  and  2c 
illustrate  the  buildup  of  the  H  atom  flux  at  lower  al¬ 
titudes.  At  550  km,  we  see  the  buildup  is  still  largest 
at  pitch  angles  near  90°,  though  by  250  km  the  pitch 
angle  distributions  for  both  protons  and  H  atoms  are 
isotropic.  Through  this  range  of  altitudes  the  models 
have  remained  in  excellent  agreement.  At  250  km  the 
MC  proton  result  is  low  compared  to  the  CSDA  and 
LT  for  the  pitch  angle  nearest  to  90°  (cosine  ==  —0.05). 
This  difference  is  apparently  from  a  slight  numerical  os¬ 
cillation  in  the  altitude  dependence  of  the  MC  results 
for  this  particular  pitch  angle.  Down  at  152  km  (Figure 


2d),  the  fluxes  remain  isotropic  and  the  models  remain 
in  good  agreement  except  near  90°  where  significant 
fall  off  of  the  flux  is  observed.  The  CSDA  shows  the 
largest  decrease  followed  by  the  LT  and  then  the  MC 
results.  By  118  km  (Figure  2e),  the  results  at  pitch 
angles  closer  to  90°  are  showing  large  decreases  in  the 
fluxes  as  well  as  some  separation  between  the  models. 
On  the  other  hand,  hearer  to  180°,  the  fluxes  are  still 
large  and  all  three  models  remain  in  good  agreement. 
As  was  seen  earlier,  when  there  is  a  significant  reduction 
in  the  flux,  it  is  the  CSDA  which  shows  the  greatest  de¬ 
crease.  At  any  particular  pitch  angle,  the  LT  and  MC 
models  remain  in  reasonable  agreement  until  the  flux 
has  decreased  about  3  orders  of  magnitude  from  what 
it  was  at  higher  altitudes.  At  110  km  (Figure  2f),  some 
separation  between  models  is  occurring  at  all  pitch  an¬ 
gles  and  by  106  km  (Figure  2g)  both  the  CSDA  and  MC 
results  are  clearly  falling  below  those  of  the  LT  model. 

In  Figures  3a  and  3b,  we  show  the  energy  integrated 
differential  flux  at  one  pitch  angle  (cosine  =  —0.95)  as 
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Cosine  of  Pilch  Angle  Cosine  of  Pitch  Angle 

Figure  2a.  The  differential  flux  integrated  over  energy  versus  the  cosine  of  the  pitch  angle  for 
(left)  protons  and  (right)  hydrogen  atoms  at  an  altitude  of  670  km.  The  error  bars  are  the  Monte 
Carlo  results,  the  long  dashes  are  the  linear  transport  results,  and  the  dotted  curves  are  the  results 
from  the  continuous  slowing-down  approximation.  The  incident  proton  flux  is  a  Maxwellian  with 
a  characteristic  energy  of  8  keV  and  a  total  incident  energy  flux  of  0.5  ergs  cm_2s_1. 


Cosine  of  Pitch  Angle  Cosine  of  Pitch  Angle 

Figure  2b.  Same  as  Figure  2a  for  an  altitude  of  550  km. 


Cosine  of  Pitch  Angle  Cosine  of  p,tch  Ang|e 

Figure  2c.  Same  as  Figure  2a  for  an  altitude  of  250  km. 
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Cosine  of  Pilch  Angle 

Figure  2g.  Same  as  Figure  2a  for  an  altitude  of  106  km. 


a  function  of  altitude.  As  before,  we  see  that  all  three 
models  are  in  excellent  agreement  until  they  reach  the 
lowest  altitudes  where  the  fluxes  are  severely  attenu¬ 
ated.  Similar  plots  for  other  pitch  angles  show  the  same 
trend  with  the  altitude  of  severe  flux  attenuation  mov¬ 
ing  higher  as  a  pitch  angle  of  90°  is  approached.  As 
seen  above,  at  all  pitch  angles  except  one  the  LT  model 
penetrates  slightly  deeper  into  the  atmosphere  than  the 
CSDA  and  MC  models.  The  one  exception  is  the  angle 
nearest  90°  (cosine  =  —0.05)  where  the  MC  penetrates 
deeper  than  the  CSDA  and  LT. 

Given  the  excellent  agreement  between  the  differen¬ 
tial  fluxes,  similar  agreement  is  expected  between  vari- 


Figure  3a.  The  differential  flux  integrated  over  energy 
versus  altitude  for  protons  (H+)  and  hydrogen  atoms 
(H)  at  a  pitch  angle  of  161.8°  (cosine  =  -0.95).  The 
error  bare  are  the  Monte  Carlo  results,  the  long  dashes 
are  the  linear  transport  results,  and  the  dotted  curves 
are  the  results  from  the  continuous  slowing-down  ap¬ 
proximation.  The  incident  proton  flux  is  a  Maxwellian 
with  a  characteristic  energy  of  8  keV  and  a  total  inci¬ 
dent  energy  flux  of  0.5  ergs  cm~8s-1. 


ous  integrals  over  the  fluxes.  In  Figure  4  we  show  the 
hemispherically  averaged  total  flux  from  all  three  mod¬ 
els  and  again  see  excellent  agreement.  Figures  5a  and 
5b  show  the  energy  deposition  rates  where  the  LT  and 
CSDA  models  agree  to  within  5%  and  both  generally 
fall  within  the  MC  errors.  Again  at  the  lowest  altitudes 
(below  the  peak  of  the  deposition)  we  see  that  the  LT 
model  penetrates  slightly  deeper  (less  than  1  km)  into 
the  atmosphere  than  does  the  CSDA  or  MC  model.  Fi¬ 
nally,  in  Figure  6,  we  show  the  eV/ion  pair  for  vari¬ 
ous  characteristic  energies  and  see  excellent  agreement 
(within  1.5%)  between  all  three  models. 

5.  Discussion  and  Summary 

We  have  shown  that  except  at  the  lowest-altitudes 
three  models  for  proton-H  atom  transport  are  in  ex¬ 
cellent  agreement.  The  differences  between  the  mod¬ 
els  are  generally  smaller  than  the  errors  that  can  arise 
from  poorly  known  cross  sections  and/or  the  errors  typi- 


Figure  3b.  Same  as  Figure  3a  with  expanded  scale  at 
lower  altitudes. 
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Figure  4.  Hemispherically  averaged  total  flux  versus 
altitude  for  protons  (H+)  and  hydrogen  atoms  (H).  The 
incident  proton  flux  is  a  Maxwellian  with  a  character¬ 
istic  energy  of  8  keV  and  a  total  incident  energy  flux  of 
0.5  ergs  cm"2s_1. 

cally  quoted  for  geophysical  observations.  For  example, 
the  measured  cross  sections  for  proton/H  atom  colli¬ 
sions  with  N2  and  O2  have  quoted  accuracies  of  typi¬ 
cally  ±20  -  30%.  For  most  collisions  involving  atomic 
oxygen,  there  are  no  measurements  and  one  resorts  to 
“estimates”  based  on  other  cross  sections.  One  excep¬ 
tion  is  for  charge  changing  collisions  for  which  measure¬ 
ments  have  been  made  with  quoted  accuracies  of  25%, 
but  at  present  there  are  still  factors  of  2  disagreement 
between  some  of  these  measurements.  Such  uncertain¬ 
ties  can  lead  to  errors  of  comparable  magnitude  in  the 


Figure  5a.  Energy  deposition  rate  versus  altitude. 
The  boxes  are  the  Monte  Carlo  results,  the  long  dashes 
are  the  linear  transport  results,  and  the  dotted  curves 
are  the  results  from  the  continuous  slowing-down  ap¬ 
proximation.  The  incident  proton  flux  is  a  Maxwellian 
with  a  characteristic  energy  of  8  keV  and  a  total  inci¬ 
dent  energy  flux  of  0.5  ergs  cm“2s_1. 


Figure  5b.  Same  as  Figure  5a  with  expanded  scale  at 
lower  altitudes. 

model  results  [Decker  et  al. ,  1995].  Likewise,  measure¬ 
ments  of  proton  fluxes  typically  involve  low  count  rates 
and  instrumental  errors  in  the  10  to  30%  range.  Optical 
measurements  of  proton  auroras  typically  quote  errors 
larger  than  10%.  Thus,  for  the  purposes  of  comparing 
to  observations  all  three  models  are  effectively  identical. 
Further,  the  quality  of  the  agreement  between  models 
gives  us  some  confidence  that  there  are  no  major  errors 
in  the  three  computer  codes. 

The  one  exception  to  the  otherwise  excellent  agree¬ 
ment  between  the  models  is  at  the  lowest  altitudes 


Figure  6.  The  quantity  weV  per  electron-ion  pair” 
versus  characteristic  energy  of  incident  protons  given 
by  Maxwellian  distributions.  The  boxes  are  the  Monte 
Carlo  results,  the  pluses  are  the  linear  transport  re¬ 
sults,  the  crosses  are  the  results  from  the  continuous 
slowing  down  approximation,  the  “E”  labels  are  results 
from  Monte  Carlo  calculations  that  include  momentum 
transfer  in  elastic  collisions  and  have  a  minimum  energy 
of  100  eV.  The  solid  curve  shows  results  from  Monte 
Carlo  calculations  that  used  the  original  cross-section 
set  of  Kozelov  and  Ivanov. 
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(see  Figures  If,  2g,  and  3b),  where  the  proton  and  H 
atom  fluxes  arc  severely  attenuated.  While  the  dif¬ 
ferences  between  the  models  are  apparent  in  all  the 
altitude-dependent  quantities  examined,  their  funda¬ 
mental  source  is  the  differences  between  the  differential 
fluxes.  In  Figures  2a-2g  it  is  evident  that  the  altitudes 
at  which  these  differences  occur  depend  on  the  pitch 
angle  of  the  flux.  However,  if  we  consider  the  fluxes  not 
as  a  function  of  altitude  but  rather  as  a  function  of  a 
pitch  angle  dependent  collision  or  scattering  depth  such 
as 

*b  f 

TP(z,  E,ii)  =  -J2  [  —»*(*').  (7) 

O  J  P 

z 

we  find  that  the  differences  between  models  occur  at 
around  the  same  collision  depth  essentially  independent 
of  the  particular  pitch  angle.  It  is  the  dependence  of 
the  collision  depth  on  the  cosine  of  the  pitch  angle  that 
causes  the  differences  between  models  to  appear  at  dif¬ 
ferent  altitudes  for  different  pitch  angles.  We  thus  find 
that  the  model  differences  arising  at  low  altitudes  oc¬ 
cur  at  large  collision  depths  where  the  collisional  mean 
free  paths  are  small  compared  to  the  scale  heights  of 
the  neutral  constituents  and  the  altitude  step  sizes  that 
are  typically  used  in  transport  calculations.  One  conse¬ 
quence  of  this  is  that  the  algorithms  designed  to  solve 
the  transport  problem  often  have  difficulties  at  these  al¬ 
titudes  and  can  be  very  sensitive  to  the  details  of  their 
numerical  implementation.  For  the  LT  model,  we  have 
found  that  it  is  particularly  sensitive  to  the  altitude 
grid  used  in  this  region.  In  our  initial  calculations  of 
the  energy  integrated  proton  differential  flux  using  a 
nonuniform  grid  of  77  altitudes,  the  LT  results  near  a 
pitch  angle  of  180°  and  at  106  km  were  over  a  factor  of 
2  larger  than  the  MC  results.  The  calculations  shown 
here  (  Figures  2g  and  3b)  used  353  altitudes,  and  the  LT 
results  are  within  20%  of  the  MC  results.  The  general 
observed  tendency  was  that  when  too  few  grid  points 
were  used  the  LT  model  penetrated  deeper  into  the  at¬ 
mosphere  than  did  the  MC  model,  and  as  the  grid  reso¬ 
lution  was  improved  the  LT  model  approached  the  MC 
model.  Fortunately,  these  differences  have  little  effect 
on  such  issues  as  where  energy  is  deposited  in  the  atmo¬ 
sphere  or  where  the  bulk  of  the  ionization  takes  place. 
Likewise,  predictions  of  the  typical  observables:  particle 
fluxes,  ion  densities,  ion  and  neutral  temperatures,  and 
optical  emissions  are  little  affected  by  these  differences. 
Thus,  for  most  aeronomical  purposes  these  differences 
deep  in  the  atmosphere  are  of  no  consequence. 

It  is  interesting  that  when  the  LT  calculations  were 
made  with  a  grid  of  77  altitudes,  the  particle  and  en¬ 
ergy  conservation  was  quite  good.  Over  most  altitudes, 
particle  and  energy  conservation  was  better  than  1%, 
and  at  the  lowest  altitudes  both  particles  and  energy 
were  conserved  to  within  a  few  percent.  When  the  cal¬ 
culations  were  made  with  353  altitudes,  the  conserva¬ 
tion  properties  were  little  changed.  For  example,  at 
some  lower  altitudes  the  particle  conservation  was  a  lit¬ 
tle  better  by  a  couple  of  percent  and  at  others  a  little 
worse  by  a  couple  of  percent.  However,  if  we  examine 


the  energy  integrated  differential  flux  near  a  pitch  an¬ 
gle  of  180°,  and  at  106  km  a  modest  change  in  particle 
conservation  is  accompained  by  over  an  80%  change  in 
the  flux.  This  can  be  understood  by  recalling  that  what 
is  normally  done  when  testing  the  particle  and  energy 
conservation  is  to  first  sum  the  total  number  of  particles 
and  total  energy  coming  into  the  system  at  the  top  of 
the  atmosphere.  One  then  calculates  at  each  altitude 
the  total  number  of  particles  and  total  energy  arriving 
at  that  altitude  along  with  what  has  been  lost  from  the 
system  at  all  higher  altitudes.  In  this  way  a  check  is 
made  whether  any  particles  or  energy  are  lost  or  cre¬ 
ated  as  a  result  of  the  numerical  methods  being  used. 
However,  once  the  lowest  altitudes  have  been  reached 
and  the  fluxes  are  orders  of  magnitude  below  their  orig¬ 
inal  values  most  of  the  particles  and  energy  have  left  the 
system.  The  number  of  particles  that  are  left  contribute 
very  little  to  the  total  conservation,  and  thus  there  can 
be  significant  errors  in  the  fluxes  yet  no  strong  indica¬ 
tion  of  this  in  the  conservation  tests.  This  implies  that 
checking  the  overall  particle  and  energy  conservation  is 
a  necessary  but  not  a  sufficient  test  of  the  accuracy  of 
a  particle  transport  model.  On  the  other  hand,  one  can 
test  particle  and  energy  conservation  at  a  particular  al¬ 
titude  relative  to  an  altitude  that  is  fairly  close.  In  that 
case,  it  is  possible  to  get  some  measure  of  the  quality  of 
the  solution  even  when  most  of  the  particles  have  left 
the  system.  For  example,  if  we  test  particle  conserva¬ 
tion  at  a  particular  altitude  relative  to  the  altitude  one 
grid  point  higher,  we  find  for  the  LT  calculation  using 
77  altitudes  that  near  a  pitch  angle  180°  between  109 
and  104  km  the  conservation  degrades  from  1%  to  26%. 
Thus  we  do  have  some  indication  that  with  this  partic¬ 
ular  altitude  grid  there  are  problems  with  the  solution 
at  the  very  low  altitudes. 

While  the  models  have  been  brought  into  good  agree¬ 
ment,  we  should  note  that  differences  can  be  seen  be¬ 
tween  the  models  as  they  are  used  in  different  stud¬ 
ies.  This  is  perhaps  most  obvious  when  Monte  Carlo 
runs  are  made  that  include  mirroring  and  beam  spread¬ 
ing.  Kozelov  [1993]  has  shown  that  these  effects  have 
an  impact  on  the  albedo,  the  charge  fractions,  and  the 
energy  deposition.  Clearly,  when  these  effects  are  im¬ 
portant  for  a  particular  study,  neither  the  CSDA  nor 
the  present  LT  model  will  be  appropriate,  and  the  MC 
model  is  needed.  In  such  a  case,  if  the  MC  is  com¬ 
putationally  too  intensive  for  the  particular  study,  one 
could  resort  to  a  table  lookup  of  a  set  of  precalculated 
MC  results.  Further  differences  can  arise  between  the 
models  due  to  differences  in  the  collisional  processes, 
Cross  sections,  and  energy  losses  that  are  used.  To  ob¬ 
tain  the  level  of  agreement  presented  in  this  paper,  we 
had  to  take  great  care  in  making  sure  that  all  these  fac¬ 
tors  were  the  same  in  all  three  models.  In  particular,  we 
found  that  the  eV/ion  pair  is  very  sensitive  to  any  such 
differences.  This  is  consistent  with  the  discussions  by 
Edgar  et  al.  [1973]  and  Strickland  et  al.  [1993]  where 
it  is  pointed  out  that  the  shape  of  the  eV /ion  pair  ver¬ 
sus  characteristic  energy  depends  on  the  behavior  of  the 
cross  sections  and  energy  losses.  For  example,  in  Figure 
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6  we  have  included  some  additional  calculations  of  the 
eV/ion  pair.  The  symbol  “E”  labels  MC  calculations 
that  used  the  cross  section  set  used  for  this  study  but 
also  included  energy  loss  due  to  elastic  collisions  and  ex¬ 
tended  the  minimum  energy  down  to  100  eV.  Here  we 
see  how  including  another  channel  for  energy  loss  (mo¬ 
mentum  transfer  in  elastic  scattering)  leaves  less  energy 
available  for  ionization  and  hence  acts  to  increase  the 
eV/ion  pair.  The  solid  curve  gives  the  results  of  MC 
calculations  of  Kozelov  and  Ivanov  [1994]  using  their 
standard  cross  section  set.  Here  the  dramatic  difference 
comes  in  large  part  from  larger  low-energy  excitation 
cross  sections  in  the  Kozelov  and  Ivanov  cross-section 
set.  Again,  the  effect  is  to  raise  the  eV/ion  pair  at 
low  energies  as  other  collisional  processes  compete  with 
ionization  to  degrade  the  energy  of  the  protons  and  H 
atoms.  The  differences  between  cross-section  sets  arise 
because  of  the  lack  of  necessary  cross-section  measure¬ 
ments  and  thus  the  need  to  make  estimates  and  ex¬ 
trapolations  in  order  to  assemble  a  complete  set.  The 
focus  of  this  paper  has  been  on  comparing  three  models 
and  not  the  issue  of  what  cross  sections  to  use  nor  how 
well  the  models  compare  to  observations.  However,  the 
success  of  the  comparisons  does  suggest  that  our  abil¬ 
ity  to  model  actual  observations  is  presently  limited  by 
uncertainties  in  cross  sections  and  the  lack  of  suitable 
observations  rather  than  our  ability  to  solve  the  equa¬ 
tions  that  describe  the  known  physics  of  proton/H  atom 
transport. 


Appendix 

The  purpose  of  this  appendix  is  to  derive  an  imple¬ 
mentation  of  the  continuous  slowing-down  approxima¬ 
tion  (CSDA)  from  linear  transport  theory.  We  begin 
with  the  linear  transport  equations  for  protons  (P)  and 
hydrogen  atoms  (H)  as  given  by  Basu  et  al  [1990] 
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where  $p(z,E,fi)  is  the  differential  flux  of  particles  (in 
units  ofcm_2s"1eV1sr“1)  of  type  0(=  P,  H)  as  a  func¬ 
tion  of  altitude  (z),  energy  (P),  and  cosine  of  the  angle 
between  the  particle  velocity  and  the  z  axis  (/z);  the  z 
axis  is  antiparallel  (parallel)  to  the  geomagnetic  field 
line  in  the  northern  (southern)  hemisphere;  and  nQ  is 
the  concentration  of  the  neutral  specie  a.  The  cross 
sections  for  the  various  collisional  processes  are  a^^(P) 
(in  units  of  cm"2),  where  j  (=fc, 10,01)  labels  the  type 
of  collision.  The  index  k  refers  only  to  excitation  and 
ionization  type  collisions  and  10  and  01  refer  to  charge 
exchange  and  stripping  collisions,  respectively.  The  to¬ 
tal  cross  section  summed  over  types  of  collisions  is  given 
by  <ja,j e{E).  For  the  differential  cross  sections,  we  make 
the  forward  scattering  approximation  to  the  angular  de¬ 
pendence  and  assume 

=  <AE') 6 (**'  -  m)  *  [&  -  wlAE')  ~  E\  (A3) 


where  Ef  is  the  incident  particle  energy,  E  is  the  final 
energy,  and  p(E)  is  the  energy  loss  associated  with 
collision  type  j. 

Our  first  step  in  developing  a  CSDA  implementation 
is  to  break  the  calculation  into  two  pieces.  On  the  one 
hand,  we  calculate  the  relative  charge  state  composition 
of  the  flux  by  neglecting  energy  loss  and  only  consider¬ 
ing  the  conversion  of  protons  into  H  atoms  and  H  atoms 
back  into  protons.  This  is  done  by  setting  W^p(E)  =  0 
and  using  (A3)  in  performing  the  fif  and  Ef  integrations 
to  obtain 
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which  are  (in  a  slightly  different  form)  (4)  and  (5)  that 
are  used  in  calculating  the  charge  state  composition  of 
the  flux.  These  equations  are  identical  to  (17)  and  (18) 
in  the  work  of  Basu  et  al  [1990]  and  were  originally  de¬ 
rived  for  use  at  high  altitudes.  While  the  particle  fluxes 
that  are  predicted  from  these  equations  are  reasonable 
only  for  high  altitudes,  it  was  found  that  the  calculated 
flux  fractions  are  accurate  at  all  altitudes. 

On  the  other  hand,  we  also  need  a  method  to  calcu¬ 
late  the  effects  of  energy  degradation  on  the  fluxes.  This 
is  done  by  assuming  the  flux  fractions  are  known  and 
calculating  an  “average”  energy  loss  for  protons  and  H 
atoms  together.  To  develop  the  needed  equation,  we 
multiply  (Al)  and  (A2)  by  E ,  integrate  over  fi  and  E , 
and  add  the  resulting  equations  to  obtain 
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Note  that  these  are  the  loss  functions  given  in  section  2 
as  (2)  and  (3).  At  this  point,  we  make  the  continuous 
slowing-down  approximation.  That  is,  we  assume  that 
the  energy  of  a  particle  is  a  continuous  function  of  al¬ 
titude,  E(z).  This  means  that  the  stochastic  nature  of 
the  collisional  degradation  of  the  protons  and  H  atoms 
is  neglected,  and  all  particles  starting  with  the  same  en¬ 
ergy  will  lose  the  same  amount  of  energy  as  they  pass 
through  the  same  range  of  altitudes.  Thus,  if  we  divide 
into  monoenergetic,  monodirectional  streams  at  the 
top  of  the  atmosphere,  they  will  maintain  their  inde¬ 
pendent  identity  as  they  penetrate  down  through  the 
atmosphere.  Thus  the  CSDA  allows  us  to  write  in 
the  following  form: 


*T(z,E,ri  =  Y,Is(z)HE-Es)6(v-fis).  (A8) 
s 


For  each  stream,  denoted  by  S,  we  can  use  (A8)  in 
evaluating  (A6)  giving 

A*s2?s/s(2:)1 


dz 


=  -'£na{z)Is(z) 

a 


La,P(Es)fp(z,  Es,  fis) 


+  Ea,B{Es)fB{z,Es,fis) 


(A9) 


We  next  return  to  (Al)  and  (A2),  integrate  over  p, 
and  E,  add  the  resulting  equations,  and  obtain  the 
statement  of  particle  conservation  that  follows  from  the 
transport  equations,  namely, 

a  0 

—  j  dEdnti$T{z,E,ti)  =  0.  (A10) 

This  gives  for  each  stream, 


d_ 

dz 


Combining  (A9)  and  (All)  gives 


i  =0. 


(All) 


dEs  1  /  \ 

dz  fj.$ ^  *w 


La,p{Es)fp(z ,  ES,fls) 


+  LCtyH{Es)fE  (z,  Es ,  fis) 


(A12) 


which  in  section  2  is  given  as  (1).  This  then  gives  us 


an  ordinary  differential  equation  to  solve  for  the  energy 
(E)  as  a  function  of  altitude  (z). 

In  the  implementation  used  in  this  paper,  we  begin 
with  a  flux  of  downgoing  protons  at  some  boundary 
altitude.  As  mentioned,  the  energy  and  pitch  angle 
distribution  of  the  protons  is  represented  by  a  set  of 
monoenergetic,  monodirectional  protons  stream.  The 
atmosphere  is  divided  into  a  series  of  horizontal  layers, 
and  (A4),  (A5),  and  (A12)  are  applied  to  each  stream. 
The  equations  are  integrated  from  altitude  to  altitude 
using  an  explicit  Euler  method  with  (A4)  and  (A5)  be¬ 
ing  solved  assuming  the  energy  of  the  stream  is  un¬ 
changed  within  a  layer  and  (A12)  being  solved  assum¬ 
ing  the  flux  fractions  arc  unchanged  within  a  layer.  At 
each  succeeding  altitude,  the  flux  fractions  and  stream 
energy  are  updated  for  use  in  integrating  down  to  the 
next  altitude.  In  this  way,  the  energy  degradation  and 
changing  composition  of  each  stream  are  calculated  as 
the  protons  and  H  atoms  penetrate  down  from  the  top 
of  the  atmosphere. 
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Abstract.  An  experimental  campaign  to  measure  diurnal  changes  in  total  electron  content 
CTEQ  over  the  wide  latitude  range  from  approximately  50°N  to  40°S  was  carried  out  from 
March  28  through  April  1 1 , 1994,  by  monitoring  the  differential  carrier  phase  from  the  U.S. 
Navy  Navigation  Satellite  System  using  a  chain  of  ground  stations  aligns  along  the 
approximate  70°W  longitude  meridian.  This  Pan-American  campaign  was  conducted  primarily 
to  study  the  day-to-day  variability  of  the  equatorial  anomaly  region.  The  experimental  plan 
included  using  die  received  values  of  TEC  from  the  chain  of  stations  to  construct  profiles  of 
electron  density  versus  latitude  using  tomographic  reconstruction  techniques  and.  then,  to 
compare  these  reconstructions  against  a  theoretical  model  of  the  low-latitude  ionosphere.  The 
diurnal  changes  in  TEC  along  this  latitude  chain  of  stations  showed  a  high  degree  of  variability 
from  day  to  day,  especially  during  a  magnetic  storm  which  occurred  near  the  beginning  of  the 
campaign.  The  equatorial  anomaly  in  TEC  showed  large  changes  in  character  in  the  two 
hemispheres,  as  well  as  differences  in  magnitude  from  day  to  day.  The  latitudinal  gradients  of 
TEC,  especially  in  the  lower  midlatitudes,  also  showed  large  differences  between  magnetic 
storm  and  quiet  conditions.  Comparisons  of  the  TEC  data  with  the  theoretical  model  illustrate 
tiie  sensitivity  of  the  model  calculations  to  changes  in  magnetic  E  x  B  drift  and  serves  to 
validate  the  strong  influence  that  these  drifts  have  on  the  formation  and  the  strength  of  the 
equatorial  anomaly  regions. 

Introduction 

The  equatorial  anomaly  region  is  characterized  by 
the  highest  values  of  peak  electron  density  and  total 
electron  content  in  the  worldwide  ionosphere.  These 
large  values  are  primarily  due  to  the  daytime  upward  E 
x  B  drift  at  (he  magnetic  equator  which  drives  the 
ionization  to  higher  altitudes.  As  the  layer  is  lifted  to 
higher  altitudes,  it  diffuses  down  magnetic  field  lint* 
that  connect  to  the  ionosphere  north  and  south  of  the 
equator,  causing  crests  at  latitudes  approximately  +/. 

15°  from  the  magnetic  equator  [Hanson  and  Moffett, 

1966,  Anderson,  1973].  This  phenomenon  has  been 
described  as  a  fountain  effect  A  discussion  on  the 
origin  of  the  terms  related  to  (he  equatorial  fountain 
effect  is  presented  in  the  appendix.  The  day-to-day 
variability  of  the  equatorial  ionosphere  can  be  very 
large,  owing  to  changes  in  the  strength  of  the  E  x  B 
drift  and  from  differences  in  the  day-to-day  strength  of 
both  the  meridional  and  zonal  component  of  neutral 
winds  [Sterling  et  al„  1969], 
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While  many  measurements  of  the  equatorial  and 
low-latitude  ionosphere  have  been  made  by  various 
techniques,  beginning  with  bottomside  ionosondes; 
then,  by  incoherent  scatter  radars  at  Jicamarca,  Peru, 
and  at  Kwajalein  Island;  later,  by  topside  sounders;  and 
finally,  by  a  dual-frequency,  dispersive  radar  on  the 
TOPEX/POSEIDON  satellite,  continuous 
measurements  of  the  ionosphere  over  a  wide  latitude 
range,  in  a  single  longitude  sector,  have  not  been 
generally  available.  Deshpande  et  al.  [1977]  first 
showed  the  day-to-day  variability  of  the  low-latitude 
total  electron  content  (TEC)  and  later,  Rastogi  and 
Klobuchar  [1990]  showed  examples  of  the  day-to-day 
variability  in  TEC  over  the  Indian  subcontinent.  More 
recently,  Whalen  [1993,  1996]  used  multiple 

ionospheric  sounders  to  illustrate  the  variability  in 
latitudinal  profiles  of  maximum  F2  region  electron 
densities  through  -t-/-30°  dip  latitude.  Model  studies  by 
Klobuchar  et  al.  [1991]  showed  that  different  E  x  B 
drifts  could  significantly  change  the  TEC  and  the 
electron  density  at  the  F2  region  peak  over  a  wide  range 
of  latitudes,  but  little  actual  diurnal  data  over  a  wide 
range  of  latitudes  have  been  available  to  check  the 
validity  of  this  model.  This  Pan-American  campaign 
was  designed  to  correct  this  deficiency. 

In  order  to  measure  the  day-to-day  variability  of  the 
ionosphere  over  a  wide  range  of  latitudes,  a  chain  of 
stations  was  set  up  along  the  approximate  70°W 
longitude  meridian  to  make  measurements  of  TEC 
using  die  dual-frequency  coherent  signals  from  the  U.S. 
Navy  Navigation  Satellite  System  (NNSS).  Figure  1 
illustrates  the  station  chain.  The  sites  were  located  at 
Hanscom  Air  Force  Base,  Massachusetts,  on  the  islands 
of  Bermuda  and  Puerto  Rico,  as  well  as  at  Merida, 
Venezuela;  La  Paz,  Bolivia,  and  Tucuman,  Argentina. 
The  planned  site  at  Tabatinga,  Brazil,  was  not  operated 
owing  to  logistical  difficulties  and  resulted  in  a  data 
gap  over  the  geomagnetic  equator.  Also  shown  in 
Figure  1  are  typical  NNSS  satellite  passes,  both  a 
northbound  and  a  southbound  one,  to  illustrate  that  the 
station  chain  was  nearly  ideal  for  such  measurements. 

Ionospheric  Measurements 

The  NNSS  satellites  transmit  coherent  radio  beacons 
at  150  and  400  MHz.  The  space  vehicles  are  placed  in 
near-circular  polar  orbits  at  altitudes  of  approximately 
1100  km  [Newton,  1967].  Evans  and  Holt  [1973] 
showed  that  these  satellites  are  ideally  suited  to  studies 
of  latitudinal  variations  of  ionospheric  behavior  by 
using  NNSS  measurements  to  observe  the  location  of 
the  midlatitude  trough  from  Millstone  Hill, 
Massachusetts.  More  recently,  tomographic 
reconstruction  studies  have  utilized  NNSS  satellite 
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measurements  in  an  effort  to  construct  profiles  of 
electron  density  versus  latitude  [Kersley  et  al.,  1993; 
Pakula  et  al.,  1995].  In  this  study,  there  were  six 
NNSS  satellites  transmitting  on  the  operational 
navigation  frequencies  during  the  approximate  2-week 
period  of  the  campaign,  making  it  possible  to  record  a 
high  elevation  pass  from  one  of  the  NNSS  satellites 
from  the  chain  of  stations  approximately  every  2  hours. 
To  determine  the  measurement  of  TEC  along  the  70° 
meridian,  the  relative  differential  carrier  pha<p 
recorded  at  all  sites  was  combined,  converted  to 
equivalent  vertical  TEC,  and  changed  into  absolute 
TEC  values  using  a  visual,  multistation  version  of  the 
two-station  method  developed  by  Leitinger  et  al 
[1975], 

The  campaign  data  are  complemented  with 
simultaneous  measurements  of  TEC  from  dual¬ 
frequency  Global  Positioning  System  (GPS)  and  from 
vertical  TEC  measured  with  the  TOPEX/POSEIDON 
satellite.  GPS  dual-frequency  data  recorded  at 
Arequipa,  Peru;  Santiago,  Chile,  and  the  island  of 
Bermuda  were  obtained  from  the  International  GPS 
Service  for  Geodynamics  network  that  is  managed  by 
the  Jet  Propulsion  Laboratory  [International 
Association  of  Geodesy,  1995].  For  this  data  set, 
absolute  measurements  of  TEC  were  mqrte  by 
combining  the  differential  phase  and  differential  group 
delay  measurements.  The  slant  TEC  values  then  were 
converted  to  equivalent  vertical  TEC  at  400  km.  In 
addition,  coincident  measurements  from  the 
TOPEX/POSEIDON  satellite  are  also  presented  in  this 
study.  The  TOPEX  satellite  mission  is  dedicated  to  the 
study  of  ocean  topography.  The  onboard  dual¬ 
frequency  (5.2  and  13.6  GHz)  altimeter  provides  near- 
global  vertical  TEC  measurements  over  the  ocean  areas 
[Intel,  1994]. 

The  measurements  were  made  near  solar  minimum 
in  late  March  and  early  April  1994.  Geomagnetic 
activity  levels  were  considered  quiet  in  the  first  few 
days  of  the  campaign.  A  major  magnetic  storm  iv»gnn 
an  April  2,  and  the  remainder  of  the  period  of 
observations  had  at  least  moderate  geomagnetic 
activity.  Data  recorded  during  the  campaign  thus 
allowed  the  opportunity  to  study  the  latihiHinql 
gradients  of  the  ionosphere  during  both 
geomagnetically  quiet  and  disturbed  conditions.  Figure 
2  illustrates  the  planetary  3-hour  Kp  indices  recorded 
during  the  campaign  period. 

Experimental  Results 

An  example  of  equivalent  vertical  TEC  from  the  full 
chain  of  six  stations  is  shown  in  Figure  3.  In  this 
figure,  the  TEC  is  plotted  from  approximately  50°N  to 
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40°S,  a  latitude  span  of  90°  encompassing  both  sides  of 
the  equatorial  region.  At  this  longitude  sector,  the 
magnetic  equator  is  located  at  approximately  11°S 
geographic  latitude.  Figure  3  shows  a  classic  case  of 
two  clearly  defined  anomaly  regions  spaced 
approximately  10  to  15°  away  from  the  geomagnetic 
equator  in  both  the  northern  and  southern  hemispheres. 
Note  that  the  two  stations  from  which  the  southern 
anomaly  peak  in  equivalent  vertical  TEC  is  inferred 
give  somewhat  different  TEC  values,  though  they  both 
show  the  peak  at  the  same  latitude.  This  difference  in 
equivalent  vertical  TEC  is  due  to  errors  in  converting 
from  the  measured  slant  TEC  to  an  equivalent  vertical 
TEC  when  looking  through  gradients  which  are  a 
function  of  elevation  angle.  This  problem  has  been 
discussed  by  Tsedilina  et  al.  [1994]. 

To  assess  the  accuracy  of  the  TEC  measurements 
maHp.  during  the  campaign,  equivalent  vertical  TEC 
from  GPS  receivers  located  at  Santiago,  Chile; 
Arequipa,  Peru,  and  the  island  of  Bermuda,  together 
with  vertical  TEC  measurements  made  via  the 
TOPEX/POSEIDON  satellite,  were  compared  with 
campaign  data.  Figure  4  illustrates  two  of  these 
comparisons.  The  GPS  measurements  shown  in  this 
figure  represent  simultaneous  measurements  of 
equivalent  vertical  TEC  at  400  km.  Since  GPS 
satellites  are  located  at  an  altitude  greater  than  20,000 
km  and  the  NNSS  satellites  are  at  approximately  1100 
km,  simultaneous  measurements  from  the  two  satellite 
systems  are  relatively  short  in  local  time.  The  TOPEX 
measurements  in  this  illustration  are  restricted  to  times 
when  the  TOPEX  ground  track  comes  within  1  hour 
local  time  of  the  campaign  measurements.  The  TOPEX 
satellite  is  located  at  approximately  1300  km  and 
provides  vertical  measurements  of  TEC  to  that  height. 
Figure  4  (left)  illustrates  close  agreement  between  the 
NNSS  and  GPS  measurements  recorded  on  day  93 
(April  3)  at  approximately  2320  IJT.  This  plot  also 
shows  that  die  TOPEX  and  NNSS  calculations  compare 
favorably  in  the  shape  and  geographic  location  of  the 
southern  anomaly  peak.  The  differences  occur  in  the 
slight  shift  in  the  anomaly  region  surrounding  the 
southern  anomaly  peak  and  the  magnitude  at  the  peak. 
The  TOPEX  calculations  illustrate  a  peak  that  is 
approximately  6  TEC  units  smaller  than  the  NNSS 
measurements.  This  difference  could  represent  an  error 
in  the  slant  to  vertical  conversion  process  used  for  the 
NNSS  data.  In  Figure  4  (right),  similar  measurements 
of  TEC  are  observed  from  all  three  observation  types 
during  an  NNSS  pass  occurring  near  1030  UT  on  day 
95  (April  5). 

The  PAN  AM  campaign  provided  the  unique 
opportunity  to  observe  the  full  diurnal  development  and 
decay  of  die  equatorial  anomaly.  Figure  5  illustrates 
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measurements  from  a  series  of  satellite  passes  for  4 
consecutive  days  during  the  campaign.  Note  that  the 
approximate  local  time  of  each  pass  is  printed  inside 
the  left  y  axis  and  that  each  pass  is  offset  by  15  TEC 
units  for  each  hour  of  local  time  after  noon  to  allow  a 
visual  separation  between  passes.  In  this  figure,  a 
classic  response  to  the  E  x  B  drift  is  portrayed,  where 
die  anomaly  begins  to  develop  by  late  morning  but  does 
not  peak  until  near  2000  hours  local  time.  Not  all  six 
stations  were  operating  for  these  days,  but  the  general 
behavior  of  the  equatorial  anomaly  and  the  differences 
in  the  shape  of  the  anomaly  among  the  4  days  can  still 
be  seen  in  Figure  5.  Note  that  the  magnitude,  width, 
and  symmetry  of  the  latitudinal  anomaly  peaks  in  the 
southern  and  northern  hemisphere  exhibit  large 
variation  over  this  4-day  period.  In  particular,  there  is 
little  evidence  of  the  northern  anomaly  peak  on  day  94 
(April  4)  at  local  noon,  while  there  is  clear  evidence  of 
it  on  days  95  through  97.  Also,  the  slope  on  the 
northern  side  of  the  anomaly  peak  at  local  noon  differs 
greatly  over  the  4  days.  Finally,  the  magnitudes  of  the 
southern  anomaly  peaks  are  much  higher  on  days  94 
and  96  than  on  days  95  and  97.  The  postmidnight  and 
early  morning  passes,  although  not  included  in  Figure 
5,  show  no  signs  of  anomaly  features.  They  are 
descriptive  of  the  downward  motion  of  the  vertical  E  x 
B  drift,  where  the  ionization  is  driven  to  lower  heights 
to  a  region  of  greater  loss.  Although  this  data  set  is 
very  informative,  the  data  gaps  near  the  magnetic 
equator  seriously  limit  our  ability  to  mnfc*  precise 
determinations  of  the  location  and  shape  of  the 
anomaly,  particularly  in  the  northern  hemisphere. 
These  data  gaps  are  due  to  the  lack  of  a  station  near  the 
magnetic  equator,  together  with  intermittent  power 
difficulties  experienced  at  the  La  Paz,  Bolivia,  station. 

The  dramatic  day-to-day  variability  illustrated  in  the 
campaign  measurements  is  consistent  with  earlier 
research.  In  particular,  Whalen  [1993,  1996] 
reconstructed  latitudinal  profiles  of  the  maximum  F 
region  electron  density  for  30  consecutive  days  in 
September  1958  by  combining  ionograms  from  a  chain 
of  ionosondes  located  in  the  western  hemisphere 
equatorial  region.  His  results  illustrate  apparent  day- 
to-day  differences  in  the  symmetry  and  magnitude  of 
the  equatorial  anomaly  peaks.  He  also  found 
asymmetry  in  the  time  of  the  enhancements,  where  the 
anomaly  peak  occurred  earlier  in  the  north  than  in  the 
south  for  most  days  in  the  observation  period. 

The  high  degree  of  variation  in  the  development  of 
flie  equatorial  anomaly  illustrated  in  Figure  5  is 
primarily  due  to  the  large  differences  in  the  vertical 
electrodynamic  lift  at  the  equator.  Drift  variations  have 
been  shown  to  have  a  large  impact  on  the  latitudinal 
location  and  amplitude  of  the  equatorial  anomaly  peaks 
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[Klobuchar  et  al.,  1991].  The  symmetry  of  both  the 
amplitude  and  position  of  the  equatorial  anomaly  peaks 
between  the  northern  and  southern  hemisphere  is 
largely  a  neutral  wind  effect.  The  next  section  focuses 
on  some  of  these  relationships  by  comparing  campaign 
data  to  a  theoretical  model  of  the  low-latitude 
ionosphere. 


Model  Comparisons 

The  data  recorded  during  the  campaign  provides  a 
unique  database  for  model  validation.  In  this  section,  a 
comparison  of  TEC  measurements  with  calculations 
from  the  Phillips  Laboratory  Global  Theoretical 
Ionospheric  Model  (GTTM)  [Anderson  et  al.,  1996],  is 
presented.  The  F  region  portion  of  this  model 
numerically  solves  the  0+  continuity  equation  to 
determine  0+  densities  as  a  function  of  altitude, 
latitude,  longitude,  and  local  time.  The  model  requires 
a  variety  of  geophysical  inputs  that  include  a  neutral 
atmosphere,  neutral  winds,  ion  and  electron 
temperatures,  and  E  x  B  drift  velocities.  The  standard 
model  inputs  and  calculations  used  in  the  low-latitude 
model  are  described  by  Preble  et  al.  [1994]  and 
Anderson  et  al.  [1996].  The  GT1M  model  is  flexible  in 
that  the  inputs  can  be  modified  to  test  the  sensitivity  of 
tiie  ionosphere  to  any  one  or  more  of  these  parameters. 

Figures  6a  through  6d  illustrate  a  comparison  of  the 
data  with  a  sequence  of  model  calculations.  Figure  6a 
shows  the  equivalent  vertical  TEC  measurements 
recorded  on  April  6.  For  clarity,  the  approximate  local 
time  is  printed  near  the  left  y  axis  and  each  curve  is 
offset  by  a  factor  of  40  TEC  units.  These  data  show  a 
well-defined  equatorial  anomaly  at  local  noon.  It  peaks 
near  1900  local  time  and  begins  to  show  signs  of  decay 
at  2300  hours.  Figures  6b-6d  represent  model 
calculations  for  comparable  geomagnetic,  seasonal,  and 
solar  conditions.  The  results  shown  in  Figure  6b  are 
based  on  a  climatological  vertical  £  x  B  drift  pattern  for 
solar  moderate  conditions  [Fejer,  1981].  Note  that  the 
local  time  development  and  the  dip  latitude  positions  of 
the  peaks  are  realistic,  while  the  magnitudes  of  the 
peaks  are  ™wh  smaller  than  those  exhibited  in  the  data 
measurements.  Initial  efforts  to  increase  the  magnitude 
of  the  equatorial  anomaly  peaks  were  made  by 
increasing  the  E  x  B  drift  The  results  of  those  efforts 
(not  shown)  illustrate  larger  peaks  but  with  the 
anomaly  pushed  out  to  higher  latitudes.  The  base  run 
of  the  model  shown  in  Figure  6b  used  a  temperature 
model  that  was  originally  developed  for  the  midlatitude 
region  [Strobel  and  McElroy,  1970].  Figure  6c 
represents  the  model  calculations  that  include  a  more 
appropriate  low-latitude  temperature  model  [Brace  and 
Theis,  1981].  This  modification  produces  sharper 
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anomaly  peaks  and  more  realistic  slopes  poleward  cf 
the  anomaly  peaks.  The  dip  latitude  positions  of  the 
anomaly  peaks  are  accurately  maintained,  but  the 
magnitude  of  the  peaks  are  still  much  smaller  than  the 
data  calculations.  The  E  x  B  vertical  drift  velocities 
used  in  the  model  are  based  on  drift  measurements 
made  at  Jicamarca,  Peru.  In  general,  they  are  applied 
in  the  GTIM  model  with  no  altitude  dependence.  Work 
by  Pingree  and  Fejer  [1987]  and  Su  et  al  [1995] 
indicate  that  the  altitude  variations  of  the  vertical  drift 
velocities  are  important  in  the  development  of  the 
equatorial  anomaly.  Figure  6d  illustrates  model  results 
when  a  simple  linear  height  variation  in  the  vertical 
drift  is  incorporated  into  the  calculations.  This 
modification  produces  more  accurate  peak  shapes  in 
both  hemispheres  and  realistic  slopes  to  higher 
latitudes  within  +/-200  dip  latitude.  The  unrealistic 
change  induced  at  latitudes  greater  than  +/-20 0 
indicates  that  the  drift  gradient  needs  to  be  refined. 
Although  the  magnitudes  of  the  anomaly  peaks  have 
been  increased,  they  are  still  lower  than  the  data 
measurements.  To  provide  a  more  quantitative 
measure  of  the  differences  shown  in  Figure  6,  plots  cf 
data  versus  model  calculations  at  1200, 1900,  and  2300 
hours  local  time  are  provided  in  Figure  7.  In  this 
figure,  it  is  evident  that  the  model  calculations  that 
include  the  linear  height  variation  in  the  drift  velocities 
best  replicate  the  features  shown  in  the  data.  This  is 
most  apparent  at  1900  hours  in  the  northern 
hemisphere  and  in  the  shape  of  anomaly  in  the  modeled 
results  at  1200  and  2300  hours. 

The  comparisons  illustrated  in  Figures  6  and  7 
illustrate  the  sensitivity  of  the  equatorial  anomaly  to  ion 
and  electron  temperatures  and  E  x  B  drift  velocities. 
Although  the  data  clearly  show  asymmetries  that  are 
likely  due  to  die  neutral  winds,  it  is  beyond  the  scope  of 
this  paper  to  replicate  this  asymmetry  with  the  model. 
Future  modeling  studies  will  include  efforts  to 
investigate  neutral  wind  effects  and  to  refine  the  E  x  B 
drift  gradient  applied  in  the  model.  In  this  study, 
GTIM  has  shown  the  capability  to  reproduce  some  of 
the  major  features  of  measured  data  in  the  low-latitude 
region,  even  under  geomagnedcally  disturbed 
conditions. 

Conclusions 

The  first  attempt  at  measuring  the  day-to-day 
variability  of  TEC  over  a  large  geographic  latitude 
range  has  been  successful.  The  database  collected 
illustrates  the  large  day-to-day  variability  in  the 
occurrence,  location,  and  amplitude  of  the  equatorial 
anomaly.  A  surprising  feature  uncovered  by  this  study 
is  that  the  TEC  values  are  so  low  in  the  latitude  range 
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greater  than  +/-400  and  that  the  latitudinal  gradients  in 
TEC  are  steeper  on  the  poleward  edge  of  the  southern 
anomaly  peak  than  on  the  poleward  edge  of  the 
northern  anomaly  peak.  The  multistation  data 
measurement  technique  has  been  validated  by 
simultaneous  measurements  of  TEC  from  dual- 
frequency  GPS  and  dual-frequency  altimeter 
measurements  from  the  TOPEX/POSEIDON  satellite. 

Model  comparisons  with  observations  have 
illustrated,  for  the  first  time,  the  sensitivity  of  the 
equatorial  anomaly  to  changes  in  ion  and  electron 
temperatures  and  to  vertical  E  x  B  drift  velocities.  In 
particular,  the  results  show  that  the  altitude  variations 
of  vertical  drift  velocities  have  a  significant  impact  on 
the  development  of  the  equatorial  anomaly.  Additional 
modeling  work  is  required  to  further  test  variations  in 
drift  gradients  and  to  investigate  neutral  wind  effects  to 
determine  how  closely  the  model  can  fit  the 
experimental  data  over  a  wide  latitude  range. 

One  of  the  prime  purposes  of  the  campaign  was  to 
develop  an  equatorial  database  to  be  used  in 
ionospheric  tomography  reconstructions.  The  study 
demonstrates  that  the  absence  of  a  station  near  the 
geomagnetic  equator  significantly  hampers  the 
reconstruction  of  the  equatorial  ionosphere.  Future 
tomographic  campaign  organizers,  take  note! 

Appendix:  Origin  of  the  Term  “equatorial 
fountain” 

Many  authors,  including  Balan  and  Bailey  [1996], 
have  shown  that  the  equatorial  fountain  greatly 
increases  both  peak  electron  density  and  TEC 
throughout  the  low-latitude  ionosphere;  but  who 
originated  this  very  apt  expression  for  the  source  of  the 
largest  ionization  in  the  world?  The  term  “fountain 
effect”  was  quoted  by  Hanson  and  Moffett  [1966, 
p.5560],  as  “Martyn  (1955)  envisages  a  ‘fountain 
effect’  in  which  the  ionization  is  lifted  upward  at  low 
latitudes  and  then  deposited  at  higher  latitudes  by 
diffusion  along  the  field  lines.”  However,  the  term 
“fountain  effect”  is  not  found  in  the  work  by  Martyn 
[1955]. 

Rush  and  Richmond  [1973,  p.1171]  stated  “Duncan 
(1960)  has  termed  this  the  ‘fountain  theory’.” 
However,  the  term  “fountain  theory”  is  not  found  in  the 
work  by  Duncan  [I960].  After  an  extensive  search,  it 
was  discovered  that  Wright  [1962]  was  the  person  who 
coined  this  term.  Wright  [1962,  p.7],  in  his  discussion 
of  the  mechanism  producing  the  equatorial  anomaly  in 
F  region  ionization,  stated  “We  propose  the  term 
‘Equatorial  Fountain’  as  a  concise  term  for  these 
processes.”  Thus  the  historical  record  is  now  set 
straight. 
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Figure  1.  Map  of  stations  along  the  •  70°  W  longitude 
meridian  used  in  the  campaign. 

Figure  2.  Planetary  3-hour  Kp  indices  recorded  during  the 
campaign  period. 
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Figure  3.  Equivalent  vertical  total  electron  content  (TEC) 
from  the  full  chain  of  six  stations  at  approximately  1600  UT 
on  April  6, 1994. 

Figure  4.  Comparison  of  TEC  measurements  from  the  U.S. 
Navy  Navigation  Satellite  System  (NNSS),  TOPEX,  and 
Global  Positioning  System  (GPS)  satellite  systems  at 
approximately  2320  hours  UT  on  April  3,  1994  (left)  and  on 
April  5  at  approximately  1030UT  (right). 

Figure  5.  Diurnal  development  of  the  equatorial  anomaly 
from  1200  to  2400  hours  local  time  for  April  4-7, 1994. 

Figure  6.  (a)  Equivalent  vertical  TEC  measurements  of  April 
6, 1994.  (b)  Model  calculations  based  on  climatological  E  x  B 
drift  and  original  temperature  model  (see  text),  (c)  Model 
calculations  with  a  more  appropriate  low-latitude  temperature 
model,  (d)  Model  calculations  with  a  more  appropriate  low- 
latitude  temperature  model  and  an  E  x  B  drift  height  gradient. 

Figure  7.  Comparison  of  data  versus  model  calculations  at 
1200, 1900,  and  2300  hours  local  time. 
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Figure  6.  (a)  Equivalent  vertical  TEC  measurements  of  April  6.  1994.  (b)  Model  calculations  based  on  climatological  E  x  B  drift  and 
original  temperature  model  (see  text),  (c)  Model  calculations  with  a  more  appropriate  low-latitude  temperature  model,  (d)  Model 
calculations  with  a  more  appropriate  low-latitude  temperature  model  and  an  E  x  B  drift  height  gradient. 

Figure  7.  Comparison  of  data  versus  model  calculations  at  1200, 1900,  and  2300  hours  local  time. 
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Geomagnetic  Latitude 
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